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Abstract. We determine the dynamical distance D, inclination i, mass-to-light ratio M/L and the intrinsic orbital 
structure of the globular cluster ui Cen, by fitting axisymmetric dynamical models to the ground-based proper 
motions of van Leeuwen et al. and line-of-sight velocities from four independent data-sets. We bring the kinematic 
measurements onto a common coordinate system, and select on cluster membership and on measurement error. 
This provides a homogeneous data-set of 2295 stars with proper motions accurate to 0.20 masyr - and 2163 stars 
with line-of-sight velocities accurate to 2 kms -1 , covering a radial range out to about half the tidal radius. 
We correct the observed velocities for perspective rotation caused by the space motion of the cluster, and show that 
the residual solid-body rotation component in the proper motions (caused by relative rotation of the photographic 
plates from which they were derived) can be taken out without any modelling other than assuming axisymmetry. 
This also provides a tight constraint on Dtani. The corrected mean velocity fields are consistent with regular 
rotation, and the velocity dispersion fields display significant deviations from isotropy. 

We model ui Cen with an axisymmetric implementation of Schwarzschild's orbit superposition method, which 
accurately fits the surface brightness distribution, makes no assumptions about the degree of velocity anisotropy 
in the cluster, and allows for radial variations in M/L. We bin the individual measurements on the plane of the sky 
to search efficiently through the parameter space of the models. Tests on an analytic model demonstrate that this 
approach is capable of measuring the cluster distance to an accuracy of about 6 per cent. Application to lo Cen 
reveals no dynamical evidence for a significant radial dependence of M/L, in harmony with the relatively long 
relaxation time of the cluster. The best-fit dynamical model has a stellar l^-band mass-to- light ratio M/Ly — 
2.5 ±0.1 M0/L0 and an inclination i = 50° ± 4°, which corresponds to an average intrinsic axial ratio of 
0.78 ±0.03. The best-fit dynamical distance D — 4.8 ±0.3 kpc (distance modulus 13.75 ±0.13 mag) is significantly 
larger than obtained by means of simple spherical or constant-anisotropy axisymmetric dynamical models, and is 
consistent with the canonical value 5.0 ± 0.2 kpc obtained by photometric methods. The total mass of the cluster 
is (2.5 ±0.3) x 10 6 Mq. 

The best-fit model is close to isotropic inside a radius of about 10 arcmin and becomes increasingly tangentially 
anisotropic in the outer region, which displays significant mean rotation. This phase-space structure may well be 
caused by the effects of the tidal field of the Milky Way. The cluster contains a separate disk-like component in 
the radial range between 1 and 3 arcmin, contributing about 4% to the total mass. 



Key words. Galaxy: globular clusters: individual: 
NGC 5139, galaxy: kinematics and dynamics 

1. Introduction 

The globular cluster u> Cen (NGC 5139) is a unique win- 
dow into astrophysics (van Leeuwen, Hughes & Piotto 
2002). It is the most massive globular cluster of our 
Galaxy, with an estimated mass between 2.4 x 1O 6 M 
(Mandushev et al. 1991) and 5.1 x 10 6 M o (Meylan et 
al. 1995). It is also one of the most flattened globular 
clusters in the Galaxy (e.g. Geyer, Nelles & Hopp 1983) 
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and it shows clear differential rotation in the line-of-sight 
(Merritt, Meylan & Mayor 1997). Furthermore, multi- 
ple stellar populations can be identified (e.g. Freeman & 
Rodgers 1975; Lee et al. 1999; Pancino et al. 2000; Bedin 
et al. 2004). Since this is unusual for a globular cluster, a 
whole range of different formation scenarios of w Cen have 
been suggested, from self-enrichment in an isolated cluster 
or in the nucleus of a tidally stripped dwarf galaxy, to a 
merger between two or more globular clusters (e.g. Icke 
& Alcaino 1988; Freeman 1993; Lee et al. 2002; Tsuchiya, 
Korchagin & Dinescu et al. 2004). 

uj Cen has a core radius of r c = 2.6 arcmin, a half- 
light (or effective) radius of Th = 4.8 arcmin and a tidal 
radius of r t — 45 arcmin (e.g. Trager, King & Djorgovski 
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1995). The resulting concentration index log(r t /r c ) ~ 1.24 
implies that ui Cen is relatively loosely bound. In com- 
bination with its relatively small heliocentric distance of 
5.0±0.2 kpc (Harris et al. 1996) 1 . This makes it is possible 
to observe individual stars over almost the entire extent 
of the cluster, including the central parts. Indeed, line-of- 
sight velocity measurements 2 have been obtained for many 
thousands of stars in the field of u> Cen (Suntzeff & Kraft 
1996, hereafter SK96; Mayor et al. 1997, hereafter M97; 
Reijns et al. 2005, hereafter Paper II; Xie, Gebhardt et 
al. in preparation, hereafter XGEA). Recently, also high- 
quality measurements of proper motions of many thou- 
sands of stars in u) Cen have become available, based 
on ground-based photographic plate observations (van 
Leeuwen et al. 2000, hereafter Paper I) and Hubble Space 
Telescope (HST) imaging (King & Anderson 2002). 

The combination of proper motions with line-of-sight 
velocity measurements allows us to obtain a dynamical 
estimate of the distance to ui Cen and study its internal 
dynamical structure. While linc-of-sight velocity observa- 
tions are in units of kms -1 , proper motions are angu- 
lar velocities and have units of (milli)arcsecyr~ 1 . A value 
for the distance is required to convert these angular ve- 
locities to kins" 1 . Once this is done, the proper motion 
and line-of-sight velocity measurements can be combined 
into a three-dimensional space velocity, which can be com- 
pared to kinematic observables that are predicted by dy- 
namical models. By varying the input parameters of these 
models, the set of model parameters (including the dis- 
tance) that provides the best-fit to the observations can 
be obtained. Similar studies for other globular clusters, 
based on comparing modest numbers of line-of-sight veloc- 
ity and proper motion measurements with simple spheri- 
cal dynamical models, were published for M3 (Cudworth 
1979), M22 (Peterson & Cudworth 1994), M4 (Peterson, 
Rees & Cudworth 1995; see also Rees 1997), and M15 
(McNamara, Harrison & Baumgardt 2004). 

A number of dynamical models which reproduce the 
line-of-sight velocity measurements for uj Cen have been 
published. As no proper motion information was included 
in these models, the distance could not be fitted and had 
to be assumed. Furthermore, all these models were lim- 
ited by the flexibility of the adopted techniques and as- 
sumed either spherical geometry (Meylan 1987, Meylan et 
al. 1995) or an isotropic velocity distribution (Merritt et 
al. 1997). Neither of these assumptions is true for co Cen 
(Geyer, Nelles & Hopp 1983; Merrifield & Kent 1990). 
Recent work, using an axisymmetric implementation of 
Schwarzschild's (1979) orbit superposition method, shows 
that it is possible to fit anisotropic dynamical models 
to (line-of-sight) kinematic observations of non-spherical 



1 Throughout the paper we use this distance of 5.0 ±0.2 kpc, 
obtained with photometric methods, as the canonical distance. 

2 Instead of the often-used term radial velocities, we adopt 
the term line-of-sight velocities, to avoid confusion with the 
decomposition of the proper motions in the plane of the sky 
into a radial and tangential component. 



galaxies (van der Marel et al. 1998; Cretton et al. 1999; 
Cappellari et al. 2002; Verolme et al. 2002; Gebhardt et 
al. 2003; Krajnovic et al. 2005). In this paper, we extend 
Schwarzschild's method in such a way that it can deal 
with a combination of proper motion and line-of-sight ve- 
locity measurements of individual stars. This allows us to 
derive an accurate dynamical distance and improve our 
understanding of the internal structure of uj Cen. 

It is possible to incorporate the discrete kinematic 
measurements of u> Cen directly in dynamical mod- 
els by using maximum likelihood techniques (Merritt & 
Saha 1993; Merritt 1993; Merritt 1997; Romanowsky & 
Kochanek 2001; Kleyna et al. 2002), but these meth- 
ods are non-linear, are not guaranteed to find the global 
best-fitting model, and are very CPU-intensive for data- 
sets consisting of several thousands of measurements. We 
therefore decided to bin the measurements instead and 
obtain the velocity moments in a set of apertures on 
the plane of the sky. While this method is (in principle) 
slightly less accurate, as some information in the data may 
be lost during the binning process, it is much faster, which 
allows us to make a thorough investigation of the param- 
eter space of uj Cen in a relatively short time. It should 
also give a good starting point for a subsequent maximum 
likelihood model using the individual measurements. 

This paper is organised as follows. In § El we describe 
the proper motion and line-of-sight velocity measure- 
ments and transform them to a common coordinate 
system. The selection of the kinematic measurements on 
cluster membership and measurement error is outlined 
in § |3 In § 0J we correct the kinematic measurements 
for perspective rotation and show that a residual solid- 
body rotation component in the proper motions can be 
taken out without any modelling other than assuming 
axisymmctry. This also provides a tight constraint on 
the inclination of the cluster. In § we describe our 
axisymmetric dynamical modelling method, and test 
it in § on an analytical model. In § we construct 
the mass model for uj Cen, bin the individual kinematic 
measurements on the plane of the sky and describe the 
construction of dynamical models that we fit to these 
observations. The resulting best-fit parameters for uj Cen 
are presented in § EI We discuss the intrinsic structure of 
the best-fit model in § [21 and draw conclusions in § ^3 

2. Observations 

We briefly describe the stellar proper motion and line-of- 
sight velocity observations of uj Cen that we use to con- 
strain our dynamical models (see Table We then align 
and transform them to a common coordinate system. 

2.1. Proper motions 

The proper motion study in Paper I is based on 100 photo- 
graphic plates of uj Cen, obtained with the Yale-Columbia 
66 cm refractor telescope. The first-epoch observations 
were taken between 1931 and 1935, for a variable star sur- 
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Table 1. Overview of the proper motions and line-of- 
sight velocity data-sets for ui Cen. The last row describes 
the four different line-of-sight velocity data-sets merged 
together, using the stars in common. The precision is es- 
timated as the median of the (asymmetric) velocity error 
distribution. If a selection on the velocity errors is applied 
(§0), the upper limit is given. For the proper motions, 
we assume a canonical distance of 5 kpc to convert from 
masyr -1 to kms -1 . 



Source Extent Observed Selected Precision 
(arcmin) (#stars) (#stars) (kms^ 1 ) 



proper motions 



Paper I 


0-30 


9847 


2295 


< 4.7 






line-of-sight velocities 






SK96 


3-23 


360 


345 


2.2 


M97 


0-22 


471 


471 


0.6 


Paper II 


0-38 


1966 


1588 


2.0 


XGEA 


0-3 


4916 


1352 


1.1 


Merged 


0-30 




2163 


< 2.0 



vey of to Cen (Martin 1938). Second-epoch plates, specif- 
ically meant for the proper motion study, were taken be- 
tween 1978 and 1983. The plates from both periods were 
compared and proper motions were measured for 9847 
stars. The observations cover a radial range of about 30 
arcmin from the cluster centre. 

2.2. Line-of-sight velocities 

We use line-of-sight velocity observations from four differ- 
ent data-sets: the first two, by SK96 and M97, from the 
literature, the third is described in the companion Paper II 
and the fourth data-set (XGEA) was kindly provided by 
Karl Gebhardt in advance of publication. 

SK96 used the ARGUS multi-object spectrograph on the 
CTIO 4 m Blanco telescope to measure, from the Ca II 
triplet range of the spectrum, the line-of-sight velocities 
of bright giant and subgiant stars in the field of 10 Cen. 
They found respectively 144 and 199 line-of-sight veloc- 
ity members, and extended the bright sample to 161 with 
measurements by Patrick Seitzer. The bright giants cover 
a radial range from 3 to 22 arcmin, whereas the subgiants 
vary in distance between 8 and 23 arcmin. From the total 
data-set of 360 stars, we remove the 6 stars without (posi- 
tive) velocity error measurement together with the 9 stars 
for which we do not have a position (see § I2.3.1J1 , leaving 
a total of 345 stars. 

M97 published 471 high-quality line-of-sight velocity 
measurements of giants in to Cen, taken with the pho- 
toelectric spectrometer CDRAVEL, mounted on the 1.5 m 
Danish telescope at Cerro La Silla. The stars in their sam- 
ple are located between 10 arcsec and 22 arcmin from the 
cluster centre. 

In Paper II, we describe the line-of-sight velocity mea- 
surements of 1966 individual stars in the field of to Cen, 



going out in radius to about 38 arcmin. Like SK96, we also 
observed with ARGUS, but used the Mg6 wavelength range. 
We use the 1589 cluster members, but exclude the single 
star for which no positive velocity error measurement is 
available. 

Finally, the data-set of XGEA contains the line-of- 
sight velocities of 4916 stars in the central 3 arcmin 
of to Cen. These measurements were obtained in three 
epochs over a time span of four years, using the Rutgers 
Imaging Fabry-Perot Spectrophotometer on the CTIO 1.5 
m telescope. During the reduction process, some slightly 
smeared out single stars were accidentally identified as two 
fainter stars. Also, contaminating light from surrounding 
stars can lead to offsets in the line-of-sight velocity mea- 
surements. To exclude (most of) these misidentifications 
(Gebhardt, priv. comm.), we select the 1352 stars with 
a measured (approximately i?-band) magnitude brighter 
than 14.5. 

2.3. Coordinate system: positions 

We constrain our dynamical models by merging all the 
above data-sets. We convert all stellar positions to the 
same projected Cartesian coordinates and align the differ- 
ent data-sets with respect to each other by matching the 
stars in common between the different data-sets. Next, we 
rotate the coordinates over the observed position angle of 
to Cen to align with its major and minor axis, and give the 
relation with the intrinsic axisymmetric coordinate system 
we assume for our models. 

2.3.1. Projected Cartesian coordinates (x",y") 

The stellar positions in Paper I are given in equatorial co- 
ordinates a and S (in units of degrees for J2000), with the 
cluster centre at a = 201?69065 and S = -47?47855. 
For objects with small apparent sizes, these equatorial co- 
ordinates can be converted to Cartesian coordinates by 
setting x" = —AacosS and y" — AS, with x" in the 
direction of West and y" in the direction of North, and 
Aa = a — ao and AS = S — So. However, this transfor- 
mation results in severe projection effects for objects that 
have a large angular diameter or are located at a large 
distance from the equatorial plane. Since both conditions 
are true for to Cen, we must project the coordinates of 
each star on the plane of the sky along the line-of-sight 
vector through the cluster centre 

x" = — rn cos <5 sin Aa, , , 

(!) 

y" = ro (sin <5 cos <5o — cos <5sin Sq cos Aa) , 

with scaling factor ro = 10800/-7T to have x" and y" in 
units of arcmin. The cluster centre is at (x" , y") = (0, 0). 

The stellar observations by SK96 are tabulated as 
a function of the projected radius to the centre only. 
However, for each star for which its ROA number (Woolley 
1966) appears in the Tables of Paper I or M97, we can re- 
construct the positions from these data-sets. In this way, 
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only nine stars are left without a position. The positions 
of the stars in the M97 data-set are given in terms of 
the projected polar radius R" in arcsec from the cluster 
centre and the projected polar angle 9" in radians from 
North to East, and can be straightforwardly converted 
into Cartesian coordinates x" and y" . For Paper II, we 
use the Leiden Identification (LID) number of each star, 
to obtain the stellar positions from Paper I. The stellar po- 
sitions in the XGEA data-set are already in the required 
Cartesian coordinates x" and y" . 

2.3.2. Alignment between data-sets 

Although for all data-sets the stellar positions are now 
in terms of the projected Cartesian coordinates (x",y"), 
(small) misalignments between the different data-sets are 
still present. These misalignments can be eliminated using 
the stars in common between the different data-sets. As 
the data-set of Paper I covers lo Cen fairly uniformly over 
much of its extent, we take their stellar positions as a 
reference frame. 

All the positions for the Paper II data-set and most 
of the positions for the SK96 data-set come directly from 
Paper I, and hence are already aligned. For the M97 and 
XGEA data-set, we use the DAOMASTER program (Stetson 
1992), to obtain the transformation (horizontal and verti- 
cal shift plus rotation) that minimises the positional dif- 
ference between the stars that are in common with those 
in Paper I: 451 for the M97 data-set and 1667 for the 
XGEA data-set. 

2.3.3. Major-minor axis coordinates {x',y') 

With all the data-sets aligned, we finally convert the stel- 
lar positions into the Cartesian coordinates (x',y'), with 
the rr'-axis and y'-axis aligned with respectively the ob- 
served major and minor axis of lo Cen. Therefore we have 
to rotate (x",y") over the position angle of the cluster. 
This angle is defined in the usual way as the angle between 
the observed major axis and North (measured counter- 
clockwise through East). 

To determine the position angle, we fit elliptic 
isophotes to the smoothed Digital Sky Survey (DSS) im- 
age of lo Cen, while keeping the centre fixed. In this way, 
we find a nearly constant position angle of 100° between 

5 and 15 arcmin from the centre of the cluster. This is 
consistent with an estimate by Seitzer (priv. comm.) from 
a J7-band image, close to the value of 96° found by White 

6 Shawl (1987), but significantly larger than the position 
angle of 91.3° measured in Paper I from star counts. 

2.3.4. Intrinsic axisymmetric coordinates (x,y,z) 

Now that we have aligned the coordinates in the plane 
of the sky (x', y') with the observed major and the major 
axis, the definition of the intrinsic coordinate system of our 
models and the relation between both becomes straight- 



forward. We assume the cluster to be axisymmetric and 
express the intrinsic properties of the model in terms of 
Cartesian coordinates (x, y, z), with the z-axis the symme- 
try axis. The relation between the intrinsic and projected 
coordinates is then given by 

x' = y, 

y' = — xcosi + zsini, (2) 
z' = — xsini — zcosi. 

The z'-axis is along the line-of-sight in the direction away 
from us 3 , and i is the inclination along which the object 
is observed, from i = 0° face-on to i = 90° edge-on. 

2.4. Coordinate system: velocities 

After the stellar positions have been transformed to a 
common coordinate system, we also convert the proper 
motions and line-of-sight velocities to the same (three- 
dimensional) Cartesian coordinate system. We centre it 
around zero (mean) velocity by subtracting the systemic 
velocity in all three directions, and relate it to the intrinsic 
axisymmetric coordinate system. 

2.4.1. Proper motions 

The proper motions (in masyr -1 ) of Paper I are given in 
the directions East and North, i.e. in the direction of — x" 
and y" respectively. After rotation over the position angle 
of 100°, we obtain the proper motion components \i x > and 
/iyi , aligned with the observed major and minor axis of 
lo Cen, and similarly, for the proper motion errors. 

2.4.2. Multiple line-of-sight velocity measurements 

In Paper II, the measured line-of-sight velocities are com- 
pared with those of SK96 and M97 for the stars in com- 
mon. A systematic offset in velocity between the differ- 
ent data-sets is clearly visible in Figure 1 of that paper. 
We measure this offset with respect to the M97 data-set, 
since it has the highest velocity precision and more than a 
hundred stars in common with the other three data-sets: 
129 with SK96, 312 with Paper II 4 and 116 with XGEA. 
As in Paper II, we apply four-sigma clipping, i.e., we ex- 
clude all stars for which the measured velocities differ by 
more than four times the combined velocity error. This 
leaves respectively 117, 284 and 109 stars in common be- 
tween M97 and the three data-sets of SK96, Paper II and 

3 In the common (mathematical) definition of a Cartesian 
coordinate system the z'-axis would point towards us, but here 
we adopt the astronomical convention to have positive line-of- 
sight away from us. 

4 In Paper II, we report only 267 stars in common with the 
data-set of M97. The reason is that there the comparison is 
based on matching ROA numbers, and since not all stars from 
Paper II have a ROA number, we find here more stars in com- 
mon by matching in position. 
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XGEA. The (weighted 5 ) mean velocity offsets of the data- 
set of M97 minus the three data-sets of SK96, Paper II and 
XGEA, are respectively -0.41 ± 0.08 kms" 1 , 1.45 ± 0.07 
kms -1 and 0.00±0.12 kms" 1 . For each of the latter three 
data-sets, we add these offsets to all observed line-of-sight 
velocities. 

Next, for each star that is present in more than one 
data-set, we combine the multiple line-of-sight velocity 
measurements. Due to non-overlapping radial coverage of 
the data-set of SK96 and XGEA, there are no stars in 
common between these two data-sets, and hence no stars 
that appear in all four data-sets. There are 138 stars with 
position in common between three data-sets and 386 stars 
in common between two data-sets. 

For the 138 stars in common between three data-sets, 
we check if the three pairwise velocity differences satisfy 
the four-sigma clipping criterion. For 6 stars, we find that 
two of the three pairs satisfy the criterion, and we select 
the two velocities that are closest to each other. For 7 
stars, we only find a single pair that satisfies the criterion, 
and we select the corresponding two velocities. Similarly, 
we find for the 386 stars in common between two data- 
sets, 13 stars for which the velocity difference does not 
satisfy the criterion, and we choose the velocity measure- 
ment with the smallest error. This means from the 524 
stars with multiple velocity measurements, for 26 stars 
(5%) one of the velocity measurements is removed as an 
outlier. This can be due to a chance combination of large 
errors, a misidentification or a binary; Mayor et al. (1996) 
estimated the global frequency of short-period binary sys- 
tems in u Cen to be 3-4%. 

As pointed out in § 2.6 of Paper II, we can use for 
the stars in common between (at least) three data-sets, 
the dispersion of the pairwise differences to calculate the 
external (instrumental) dispersion for each of the data- 
sets. In this way, we found in Paper II that the errors 
tabulated in SK96 are under- estimated by about 40% and 
hence increased them by this amount, whereas those in 
M97 are well-calibrated. Unfortunately, there are too few 
stars in common with the XGEA data-set for a similar 
(statistically reliable) external error estimate. 

In the final sample, we have 125 stars with the 
weighted mean of three velocity measurements and 373 
stars with the weighted mean of two velocity measure- 
ments. Together with the 2596 single velocity measure- 
ments, this gives a total of 3094 cluster stars with line-of- 
sight velocities. 

2.4.3. Systemic velocities 

To centre the Cartesian velocity system around zero mean 
velocity, we subtract from both the proper motion data- 
sets and the merged line-of-sight data-set the (remaining) 
systemic velocities. In combination with the cluster proper 

5 To calculate the mean and dispersion of a sample, we use 
the weighted estimators and corresponding uncertainties as de- 
scribed in Appendix A of Paper II. 



motion values from Table 4 of Paper I, we find the follow- 
ing systemic velocities 

f/j; = 3.88 ±0.41 masyr" 1 , 

f/ y y , B = -4.44 ±0.41 masyr" 1 , (3) 
v s J s = 232.02 ±0.03 kms" 1 . 

2.4.4. Intrinsic axisymmetric coordinate system 

In our models, we calculate the velocities in units of 
kms -1 . If we assume a distance D (in units of kpc), the 
conversion of the proper motions in units of masyr -1 into 
units of kms -1 is given by 

?V = 4.74 D/v and v y > = A.7AD^. (4) 

The relation between observed (u x < , v y > , v x > ) and intrinsic 
(v x , v y ,v z ) velocities is the same as in equation (JSJ, with 
the coordinates replaced by the corresponding velocities. 

In addition to Cartesian coordinates, we also describe 
the intrinsic properties of our axisymmetric models in 
terms of the usual cylindrical coordinates (R,<fi,z), with 
x = R cos 4> and y = R sin <f>. In these coordinates the 
relation between the observed and intrinsic velocities is 

v x ' = Vfjsin0 + W0COS0, 

v y i — (— w^cos0 + v^smifi) cosi + v z sini, (5) 
v z > = (—vr cos (p + sin tfi) sini + v z cos i. 

3. Selection 

We discuss the selection of the cluster members from the 
different data-sets, as well as some further removal of stars 
that cause systematic deviations in the kinematics. 

3.1. Proper motions 

In Paper I, a membership probability was assigned to each 
star. We use the stars for which we also have line-of-sight 
velocity measurements to investigate the membership de- 
termination. Furthermore, in Paper I the image of each 
star was inspected and classified according to its separa- 
tion from other stars. We study the effect of the distur- 
bance by a neighbouring star on the kinematics. Finally, 
after selection of the undisturbed cluster members, we ex- 
clude the stars with relatively large uncertainties in their 
proper motion measurements, which cause a systematic 
overestimation of the mean proper motion dispersion. 

3.1.1. Membership determination 

The membership probability in Paper I was assigned to 
each star in the field by assuming that the distribution 
of stellar velocities is Gaussian. In most studies, this is 
done by adopting one common distribution for the en- 
tire cluster. However, this does not take into account that 
the internal dispersion, as well as the relative number of 
cluster stars, decreases with radius. To better incorporate 
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disturbance by neighboring stars 
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Fig. 1. Velocity dispersion profiles, calculated along con- 
centric rings, from the proper motions of Paper I. The dis- 
persion profiles from the proper motions in the rr'-direction 
(y'-direction) are shown in the top (bottom). The error bar 
at the bottom-left indicates the typical uncertainty The 
red curves are the dispersion profiles for all 7899 clus- 
ter stars with proper motion measurements. The other 
coloured curves show how the dispersion decreases signif- 
icantly, especially in the crowded centre of u> Cen, when 
sequentially stars of class 4 (severely disturbed) to class 1 
(slightly disturbed) are removed. We select the 4415 undis- 
turbed stars of class 0. 
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Fig. 2. Proper motion dispersion profiles as in Figure ^ 
Starting with all undisturbed (class 0) cluster stars (red 
solid curve), sequentially a smaller number of stars is se- 
lected by setting a tighter limit on the allowed error in 
their proper motion measurements. The dispersion de- 
creases if the stars with uncertain proper motion measure- 
ments are excluded. This effect is significant and larger 
than the dispersion broadening due to the individual ve- 
locity errors, indicated by the red dotted curve. We select 
the 2295 stars with proper motion error smaller than 0.20 
masyr" 1 , since below this limit the kinematics stay simi- 
lar. 



these two effects, the membership probability in Paper I 
was calculated along concentric rings. 

By matching the identification numbers and the posi- 
tions of stars, we find that there are 3762 stars for which 
both proper motions and line-of-sight velocities are mea- 
sured. This allows us to investigate the quality of the mem- 
bership probability assigned in Paper I, as the separation 
of cluster and field stars is very clean in line-of-sight ve- 
locities (see e.g. Paper II, Figure 4). 

From the line-of-sight velocities, we find that of the 
3762 matched stars, 3385 are cluster members. Indeed, 
most of these cluster stars, 3204 (95%), have a member- 
ship probability based on their proper motions of at least 
68 per cent. Based on the latter criterion, the remain- 
ing 181 (5%) cluster stars are wrongly classified as field 
stars in Paper I. From the 3762 matched stars, 377 stars 
are field stars from the line-of-sight velocity data-set of 
Paper II. Based on a membership probability of 68 per 



cent, 54 (14%) of these field stars are wrongly classified 
as cluster members in Paper I. This fraction of field stars 
misclassified as cluster stars is an upper limit, since the 
obvious field stars are already removed from the proper 
motion data-set of Paper I. 

Wrongly classifying cluster stars as field stars is rela- 
tively harmless for our purpose, since it only reduces the 
total cluster data-set. However, classifying field stars as 
members of the cluster introduces stars from a different 
population with different (kinematical) properties. With 
a membership probability of 99.7 per cent the fraction 
of field stars misclassified as cluster stars reduces to 5%. 
However, at the same time we expect to miss almost 30% 
of the cluster stars as they are wrongly classified as field 
stars. Taking also into account that the additional selec- 
tions on disturbance by neighbouring stars and velocity 
error below remove (part of) the field stars misclassified 
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as cluster stars, we consider stars with a membership prob- 
ability of at least 68 per cent as cluster members. 

While for the 3762 matched stars, the line-of-sight ve- 
locities confirm 3385 stars as cluster members, from the 
remaining 6084 (unmatched) stars of Paper I, 4597 stars 
have a proper motion membership probability of at least 
68 per cent. From the resulting proper motion distribu- 
tion, we remove 83 outliers with proper motions five times 
the standard deviation away from the mean, leaving a to- 
tal of 7899 cluster stars. 

3.1.2. Disturbance by neighbouring stars 

In Paper I, each star was classified according to its sep- 
aration from other stars on a scale from to 4, from 
completely free to badly disturbed by a neighbouring 
star. In Figure ^ we show the effect of the disturbance 
on the proper motion dispersion. The (smoothed) pro- 
files are constructed by calculating the mean proper mo- 
tion dispersion of the stars binned in concentric rings, 
taking the individual measurement errors into account 
(Appendix^). The proper motions in the x'-direction give 
rise to the velocity dispersion profiles a x i in the upper 
panel. The proper motions in the ^/'-direction yield the 
velocity dispersion profiles a y i in the bottom panel. The 
red curves are the velocity dispersion profiles for all 7899 
cluster stars with proper motion measurements. The other 
coloured curves show how, especially in the crowded cen- 
tre of to Cen, the dispersion decreases significantly when 
sequentially stars of class 4 (severely disturbed) to class 1 
(slightly disturbed) are removed. We select the 4415 undis- 
turbed stars of class 0. 

The membership determination is cleaner for undis- 
turbed stars, so that above fraction of 5% of the cluster 
stars misclassified as field stars becomes smaller than 3% 
if only stars of class are selected. The velocity dispersion 
profiles a x i and ay in Figure Q ai 'e systematically offset 
with respect to each other, demonstrating that the veloc- 
ity distribution in oj Cen is anisotropic. We discuss this 
further in § EH and §E3 

3.1.3. Selection on proper motion error 

After selection of the cluster members that are not dis- 
turbed by neighbouring stars, it is likely that the sample of 
4415 stars still includes (remaining) interlopers and stars 
with uncertain proper motion measurements, which can 
lead to systematic deviations in the kinematics. Figure 
shows that the proper motion dispersion profiles decrease 
if we sequentially select a smaller number of stars by set- 
ting a tighter limit on the allowed error in their proper 
motion measurements. 

Since the proper motion errors are larger for the fainter 
stars (see also Figure 11 of Paper I), a similar effect hap- 
pens if we select on magnitude instead. The decrease in 
dispersion is most prominent at larger radii as the above 
selection on disturbance by a neighbouring star already 



removed the uncertain proper motion measurements in 
the crowded centre of u> Cen. All dispersion profiles in the 
above are corrected for the broadening due to the individ- 
ual proper motion errors (cf. Appendix EJ. The effect of 
this broadening, indicated by the dotted curve, is less than 
the decrease in the dispersion profiles due to the selection 
on proper motion error. 

Since the kinematics do not change anymore signifi- 
cantly for a limit on the proper motion errors lower than 
0.20 masyr" 1 , we select the 2295 stars with proper mo- 
tion errors below this limit. The preliminary HST proper 
motions of King & Anderson (2002) in the centre of a; Cen 
(i?' ~ 1 arcmin) give rise to mean proper motion disper- 
sion ov = 0.81 ± 0.08 masyr" 1 and oy = 0.77 ± 0.08 
masyr -1 , depending on the magnitude cut-off. In their 
outer calibration field (B! ~ 14 arcmin), the average dis- 
persion is about 0.41±0.03 masyr -1 . These values are con- 
sistent with the mean proper motion dispersion of the 2295 
selected stars at those radii (light blue curves in Figure^. 
We are therefore confident that the proper motion kine- 
matics have converged. 

The spatial distribution of the selected stars is shown 
in the top panel of Figure 0] In the two upper panels 
of Figure the distributions of the two proper motion 
components (left panels) and the corresponding errors 
(right panels) of the N sc \ = 2295 selected stars are shown 
as shaded histograms, on top of the histograms of the 
JVmem = 7899 cluster members. The selection removes the 
extended tails, making the distribution narrower with an 
approximately Gaussian shape. 

3.2. Line-of-sight velocities 

For each of the four different line-of-sight velocity data- 
sets separately, the velocity dispersion profiles of the se- 
lected (cluster) stars (§ 12.21 and Table are shown as 
solid coloured curves in Figure 01 The dotted blue curve 
is the dispersion profile of all the 4916 stars observed by 
XGEA, whereas the solid blue curve is based on the 1352 
selected stars with a measured magnitude brighter than 
14.5, showing that fainter misidentified stars lead to an 
under-estimation of the line-of-sight velocity dispersion. 
Although the dispersion profile of the M97 data-set (yel- 
low curve) seems to be systematically higher than those of 
the other data-sets, it is based on a relatively small number 
of stars, similar to the SK96 data-set, and the differences 
are still within the expected uncertainties indicated by the 
error bar. 

The solid black curve is the dispersion profile of the 
3094 stars after merging the four line-of-sight velocity 
data-sets (S I2.4.2JI . Due to uncertainties in the line-of-sight 
velocity measurements of especially the fainter stars, the 
latter dispersion profile is (slightly) under-estimated in the 
outer parts. By sequentially lowering the limit on the line- 
of-sight velocity errors, we find that below 2.0 kins" 1 the 
velocity dispersion (dotted black curve) converges. Hence, 
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Fig. 3. Velocity dispersion profiles, calculated along con- 
centric rings, for the four different line-of-sight velocity 
data-sets separately and after they have been merged. 
The blue dotted curve shows the under-estimated disper- 
sion for the XGEA data-set if also the faint stars arc in- 
cluded. From the merged data-set of 3094 stars we select 
the 2163 stars with line-of-sight velocity errors smaller 
than 2.0 kms -1 , resulting in a dispersion profile (black 
dotted curve) that is not under- estimated due to uncer- 
tain line-of-sight velocity measurements. 

we select the 2163 stars with line-of-sight velocity errors 
smaller than 2.0 kms -1 . 

The spatial distribution of the selected stars is shown 
in the bottom panel of Figure 0] In the bottom panels 
of Figure the distribution of the line-of-sight velocities 
(left) and corresponding errors (right) of the N sc \ = 2163 
selected stars are shown as filled histograms, on top of the 
histograms of the N mem — 3094 cluster members in the 
merged data-set. 

4. Kinematics 

We compute the mean velocity fields for the selected stars 
and correct the kinematic data for perspective rotation 
and for residual solid-body rotation in the proper motions. 
At the same time, we place a tight constraint on the incli- 
nation. Finally, we calculate the mean velocity dispersion 
profiles from the corrected kinematic data. 
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Fig. 4. The stars in u> Cen with proper motion measure- 
ments (top) and line-of-sight velocity measurements (bot- 
tom), that are used in our analysis. The stellar positions 
are plotted as a function of the projected Cartesian coordi- 
nates x' and y 1 , with the x'-axis aligned with the observed 
major axis and the y'-axis aligned with the observed mi- 
nor axis of qj Cen. The excess of stars with line-of-sight 
velocities inside the central 3 arcmin in the bottom panel 
is due to the XGEA data-set. 



4.1. Smoothed mean velocity fields 

The left-most panels of FigureH3show the smoothed mean 
velocity fields for the 2295 selected stars with proper 
motion measurements and the 2163 selected stars with 
line-of-sight velocity measurements. This adaptive kernel 
smoothening is done by selecting for each star its 200 near- 
est neighbours on the plane of the sky, and then calculat- 



ing the mean velocity (and higher order velocity moments) 
from the individual velocity measurements ( Appendix |A"|) . 
The contribution of each neighbour is weighted with its 
distance to the star, using a Gaussian distribution with 
zero mean and the mean distance of the 200 nearest neigh- 
bours as the dispersion. 
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Fig. 5. Histograms of measured velocities (left pan- 
els) and corresponding velocity errors (right panels). The 
proper motion components [j, x > (upper panels) and n v t 
(middle panels), in the direction of the observed major 
and minor axis of u> Cen respectively, come from the pho- 
tographic plate observations in Paper I. The line-of-sight 
velocities (lower panels) are taken from four different data- 
sets (§ I2.2|l . The shaded histograms for the iV sc i selected 
stars (§ 0J are overlayed on the histograms of the N mem 
cluster member stars. 

The upper-left panel shows the mean proper motion 
(in masyr -1 ) in the major axis ^'-direction, i.e., the hori- 
zontal component of the streaming motion on the plane of 
the sky. The colour coding is such that red (blue) means 
that the stars are moving on average to the right (left) and 
green shows the region where the horizontal component of 
the mean proper motion vanishes. Similarly, the middle- 
left panel shows the mean proper motion in the minor axis 
y'-direction, i.e. the vertical component of the streaming 
motion on the plane of the sky, with red (blue) indicating 



average proper motion upwards (downwards). Finally, the 
lower-left panel shows the mean velocity (in kms -1 ) along 
the line-of-sight z'-axis, where red (blue) means that the 
stars are on average receding (approaching) and green in- 
dicates the zero-velocity curve, which is the rotation axis 
of uj Cen. 

Apart from a twist in the (green) zero-velocity curve, 
the latter line-of-sight velocity field is as expected for 
a (nearly) axisymmetric stellar system. However, both 
proper motion fields show a complex structure, with an 
apparently dynamically decoupled inner part, far from ax- 
isymmetric. We now show that it is, in fact, possible to 
bring these different observations into concordance. 

4.2. Perspective rotation 

The non-axisymmctric features in the observed smoothed 
mean velocity fields in the left-most panels of Figure 
might be (partly) caused by perspective rotation. Because 
u> Cen has a large extent on the plane of the sky (with a 
diameter about twice that of the full moon), its substan- 
tial systemic (or space) motion (eq. produces a non- 
negligible amount of apparent rotation: the projection of 
the space motion onto the principal axis (&', y' , z') is dif- 
ferent at different positions on the plane of the sky (Feast 
et al. 1961). We expand this perspective rotation in terms 
of the reciprocal of the distance D. Ignoring the negligi- 
ble terms of order 1/D 2 or smaller, we find the following 
additional velocities 

= -6.1363 xlO- 5 x'v s z r/D masyr" 1 , 
fjP 1 ; = -6.1363 xlO- 5 y'v s z r/D masyr -1 , (6) 

v p J = 1.3790 x 1(T 3 (V//J s + y'^ s ) D kms" 1 , 

with x 1 and y' in units of arcmin and D in kpc. For the 
canonical distance of 5 kpc, the systemic motion for to Cen 
as given in eq. and the data typically extending to 20 
arcmin from the cluster centre, we find that the maxi- 
mum amplitude of the perspective rotation for the proper 
motions is about 0.06 masyr -1 and for the line-of-sight 
velocity about 0.8 kms -1 . These values are a significant 
fraction of the observed mean velocities (left panels of 
Figure [BJ) and of the same order as the uncertainties in 
the extracted kinematics (see Appendix |A"|) . Therefore, the 
perspective rotation as shown in the second column panels 
of Figure® cannot be ignored and we correct the observed 
stellar velocities by subtracting it. Since we use the more 
recent and improved values for the systemic proper motion 
from Paper I, our correction for perspective rotation is dif- 
ferent from that of Merritt et al. (1997). The amplitude 
of the correction is, however, too small to explain all of 
the complex structure in the proper motion fields and we 
have to look for an additional cause of non-axisymmetry. 

4.3. Residual solid-body rotation 

Van Leeuwen & Le Poole (2002) already showed that a 
possible residual solid-body rotation component in the 
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Fig. 6. The mean velocity fields of uj Cen corrected for perspective and solid-body rotation. The individual measure- 
ments are smoothed using adaptive kernel smoothcning. From top to bottom: The mean ground-based proper motion 
in the major axis aZ-direction and in the minor axis y'-dircction, and the mean linc-of-sight velocity. From left to right: 
Observed velocity fields of u) Cen, contribution from perspective rotation, contribution from solid-body rotation and 
the velocity fields after correcting for both. The perspective rotation is caused by the space motion of oj Cen. The 
solid-body rotation in the proper motions is due to relative rotation of the first and second epoch photographic plates 
by an amount of 0.029 masyr -1 arcmin -1 (§ I4.4|l . 



ground-based proper motions of Paper I can have an im- 
portant effect on the kinematics. The astrometric reduc- 
tion process to measure proper motions removes the abil- 
ity to observe an overall rotation on the plane of the 
sky (e.g. Vasilevskis et al. 1979). This solid-body rota- 
tion results in a transverse proper motion v t = ili?', 
with Q the amount of solid-body rotation (in units of 
masyr -1 arcmin -1 ) and R' the distance from the clus- 
ter centre in the plane of the sky (in units of arcmin). 
Decomposition of Vt along the observed major and minor 
axis yields 

fJ$" = +Vy' masyr -1 , 

ffi = -Qx' masyr -1 . [ ) 



Any other reference point than the cluster centre results 
in a constant offset in the proper motions, and is removed 
by setting the systemic proper motions to zero. Also an 
overall expansion (or contraction) cannot be determined 
from the measured proper motions, and results in a ra- 
dial proper motion in the plane of the sky. Although both 
the amount of overall rotation and expansion are in prin- 
ciple free parameters, they can be constrained from the 
link between the measured (differential) proper motions 
to an absolute proper motion system, such as defined by 
the Hipparcos and Tycho-2 catalogues (Perryman et al. 
1997; H0g et al. 2000). In Paper I, using the 56 stars in 
common with these two catalogues, the allowed amount of 
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residual solid-body rotation was determined to be no more 
than SI = 0.02±0.02 masyr -1 axcmin -1 and no significant 
expansion was found. 

As the amplitude of the allowed residual solid-body 
rotation is of the order of the uncertainties in the mean 
proper motions already close to the centre, and can 
increase beyond the maximum amplitude of the mean 
proper motions in the outer parts, correcting for it has 
a very important effect on the proper motions. We use a 
general relation for axisymmetric objects to constrain S7, 
and at the same find a constraint on the inclination. 

4.4. The amount of residual solid-body rotation 
directly from the mean velocities 

For any axisymmetric system, there is, at each position 
(x',y') on the plane of the sky, a simple relation between 
the mean proper motion in the y'-direction (/%') and the 
mean line-of-sight velocity (v z i) (see e.g. Appendix A of 
Evans & de Zeeuw 1994, hereafter EZ94). Using rela- 
tion (JSJ, with for an axisymmetric system (vr) = (v z ) = 0, 
we see that, while the mean velocity component in the in- 
direction includes the spatial term cos <j), those in the y'- 
direction and line-of-sight z'-direction both contain sin <f>. 
The latter implies that, by integrating along the line-of- 
sight to obtain the observed mean velocities, the expres- 
sions for {vyi) and (v z r) only differ by the cosi and sini 
terms. Going from (ty) to via equation we thus 
find the following general relation for axisymmetric ob- 
jects 

(v z ,)(x',y') = 4.74£ta ni {fi y ,)(x', y'), (8) 

with distance D (in kpc) and inclination i. 

This relation implies that, at each position on the plane 
of the sky, the only difference between the mean short- 
axis proper motion field and the mean line-of-sight ve- 
locity field should be a constant scaling factor equal to 
4.74 Dt&ni. Comparing the left-most middle and bottom 
panel in Figure (V^bserved)) this is far from what we see, 
except perhaps for the inner part. We ascribe this discrep- 
ancy to the residual solid-body rotation, which causes a 
perturbation of (/jL y >) increasing with x' as given in eq. JJJ. 
In this way, we can objectively quantify the amount of 
solid body rotation S7 needed to satisfy the above relation 
(jHJ , and at the same time find the best-fit value for D tan i. 

To compute uncorrelated values (and corresponding 
errors) for the mean short-axis proper motion (ny>} and 
mean line-of-sight velocity {v z >) at the same positions on 
the plane of the sky, we bin the stars with proper motion 
and line-of-sight velocity measurements in the same polar 
grid of apertures (see also Appendix |A"|) . We plot the re- 
sulting values for (v z i) against {n y >) and fit a line (through 
the origin) by minimising the \ 2 , taking into account the 
errors in both directions (§ 15.3 of Press et al. 1992). 

By varying the amount of solid-body rotation S7 and 
the slope of the line, which proportional to D tan i (eq. |HJ) , 
we obtain the A\ 2 = X 2 — Xmin contours in the left 



panel of Figure [7| The inner three contours are drawn 
at the levels containing 68.3%, 95.4% and 99.7% (thick 
contour) of a A% 2 -distribution with two degrees of free- 
dom. 6 Subsequent contours correspond to a factor of two 
increase in A\ 2 . The overall minimum Xmim indicated by 
a cross, implies (at the 68.3%-level) a best-fit value of SI = 
0.029 ± 0.004 masyr -1 arcmin -1 . This is fully consistent 
with the upper limit of SI = 0.02±0.02 masyr -1 arcmin -1 
from Paper I. 

The middle panel of Figure shows that without any 
correction for residual solid-body rotation, the values for 
(v z ') and (n y >) are scattered (open circles), while they 
are nicely correlated after correction with SI = 0.029 
masyr -1 arcmin _1 (filled circles). The resulting solid- 
body rotation, shown in the third column of Figure El 
removes the cylindrical rotation that is visible in the 
outer parts of the observed proper motion fields (first col- 
umn). After subtracting this residual solid-body rotation, 
together with the perspective rotation (second column), 
the complex structures disappear, resulting in (nearly) 
axisymmetric mean velocity fields in the last column. 
Although the remaining non-axisymmetric features, such 
as the twist of the (green) zero- velocity curve, might indi- 
cate deviations from true axisymmetry, they can also be 
(partly) artifacts of the smoothening, which, especially in 
the less dense outer parts, is sensitive to the distribution 
of stars on the plane of the sky. 

This shows that the application of eq. JHJ to the com- 
bination of proper motion and line-of-sight measurements 
provides a powerful new tool to determine the amount of 
solid body rotation. At the same time, it also allows us to 
place a constraint on the inclination. 

4.5. Constraint on the inclination 

From the left panel of Figure we obtain (at the 68.3%- 
level) a best-fit value for Z?tani of 5.6 (+1.9/-1.0) kpc. 
The right panel shows the resulting relation (solid line) 
between the distance D and the inclination i, where the 
dashed lines bracket the 68.3%-level uncertainty. If we as- 
sume the canonical value D = 5.0 ± 0.2 kpc, then the 
inclination is constrained to i = 48 (+9/-7) degrees. 

Although we apply the same polar grid to the proper 
motions and line-of-sight velocities, the apertures contain 
different (numbers of) stars. To test that this does not 
significantly influence the computed average kinematics 
and hence the above results, we repeated the analysis but 
now only include the 718 stars for which both the proper 
motions and line-of-sight velocity are measured. The re- 



For a Gaussian distribution with dispersion a, these per- 
centages correspond to the la, 2a and 3a confidence intervals 
respectively. For the (asymmetric) x 2 -distribution there is in 
general no simple relation between dispersion and confidence 
intervals. Nevertheless, the 68.3%, 95.4% and 99.7% levels of 
the ^-distribution are often referred to as the la, 2a and 3<r 
levels. 
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Fig. 7. Constraints on the amount of residual solid-body body rotation and via Dtani, on the distance D (in 
kpc) and inclination i, using the general relation JSJ for axisymmetric objects. The left panel shows the contour map 
of the goodness-of-fit parameter A\ 2 . The inner three contours are drawn at the 68.3%, 95.4% and 99.7% (thick 
contour) levels of a Ax 2 -distribution with two degrees of freedom. Subsequent contours correspond to a factor of two 
increase in A% 2 . The overall minimum is indicated by the cross. The middle panel shows the mean line-of-sight velocity 
(v z >) and mean short-axis proper motion (fj, y >) within the same polar apertures, before (open circles) and after (filled 
circles) correction for residual solid-body rotation with the best-fit value of f2 = 0.029 ± 0.004 masyr -1 arcmin -1 . The 
best-fit value for Dtani = 5.6 (+1.9/-1.0) kpc gives rise to the relation in the right panel (sold line), bracketed (at 
the 68.3%-level) by the dashed lines. Given the canonical distance of D = 5.0 ± 0.2 kpc, the dotted lines indicate the 
constraint on inclination of i = 48 (+9/-7) degrees. 



suits are equivalent, but less tightly constrained due to 
the smaller number of apertures. 

Van Leeuwen & Le Poole (2002) compared, for differ- 
ent values for the amount of residual solid-body rotation 
fl, the shape of the radial profile of the mean transverse 
component of proper motions from Paper I, with that of 
the mean line-of-sight velocities calculated by Merritt et 
al. (1997) from the line-of-sight velocity data-set of M97. 
They found that il ~ 0.032 masyr -1 arcmin -1 provides a 
plausible agreement. Next, assuming a distribution for the 
density and the rotational velocities in the cluster, they 
computed projected transverse proper motion and line- 
of-sight velocity profiles, and by comparing them to the 
observed profiles, they derived a range for the inclination i 
from 40 to 60 degrees. Their results are consistent with our 
best-fit values = 0.029 ± 0.004 masyr -1 arcmin -1 and 
i = 48 (+9/-7) degrees. Our method is based on a general 
relation for axisymmetric objects, without any further as- 
sumptions about the underlying density and velocity dis- 
tribution. Moreover, instead of comparing shapes of mean 
velocity profiles, we actually fit the (two-dimensional) 
mean velocity fields. 

In the above analysis, we assume that all of the solid- 
body rotation in the proper motion is the result of a 
(non-physical) residual from the photographic plate re- 
duction in Paper I. This raises the question what happens 
if a (physical) solid-body rotation component is present in 
u> Cen. Such a solid-body rotation component is expected 
to be aligned with the intrinsic rotation axis, inclined at 
about 48°, and therefore also present in the line-of-sight 
velocities. Except for the perspective rotation correction, 



we leave the mean line-of-sight velocities in the above anal- 
ysis unchanged, so that any such solid-body rotation com- 
ponent should also remain in the proper motion. 

Still, since we are fitting the residual solid-body ro- 
tation and the slope D tan i simultaneously, we show 
next that these parameters can become (partly) degener- 
ate. Combining eq. (JJJ) with . we obtain the best-fit 
values for D tan i and Q by minimising 

2 
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where (v°,)j and are respectively the observed 

mean line-of-sight velocity and the observed mean proper 
motion in the y'-direction, measured in aperture j of a 
total of n apertures, with their centres at x'y A(v°f 3S ) J - 
and A(fi°]' s )j are the corresponding measurement errors. 
Suppose now the extreme case that all of the observed 
mean motion is due to solid-body rotation: an amount of 
f^o residual solid-body rotation in the plane of the sky, and 
an amount of loq intrinsic solid-body rotation, around the 
intrinsic z-axis in u> Cen, which we assume to be inclined 
at io degrees. At a distance Dq, the combination yields 
per aperture (v°}' s )j = 4.74DqWo sin iox'j and ([J,ff s )j = 
(wocosio — f2o)x^. Substitution of these quantities in the 
above eq. ©, and ignoring the (small) variations in the 
measurements errors, yields that x 2 = if 

-l 



D tan i = Dq tan i$ 



1 



(10) 



UlQ COS lQ 

This implies a degeneracy between Dtani and fi, which 
left panel of Figure [7| would result in the same minimum 
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Fig. 8. Mean velocity dispersion profiles calculated along 
concentric rings. Assuming the canonical distance of 5 kpc, 
the profiles of the radial aw (green) and tangential ag* 
(red) components of the proper motion dispersion are con- 
verted into the same units of kms -1 as the profile of the 
line-of-sight velocity dispersion a z > (blue). The black hor- 
izontal lines indicate the corresponding scale in masyr -1 . 
The mean velocity error per ring is indicated below the 
profiles by the crosses. The green and red crosses mostly 
overlap, as the errors of the radial and tangential compo- 
nents are nearly similar. 

all along a curve. However, in the case the motion in 10 Cen 
consists of more than only solid-body rotation, this degen- 
eracy breaks down and we expect to find a unique mini- 
mum. The latter seems to be the case here, and we con- 
clude that the degeneracy and hence the intrinsic solid- 
body rotation are not dominant, if present at all. 

4.6. Mean velocity dispersion profiles 

In Figure we show the mean velocity dispersion profiles 
of the radial <tr/ (green) and tangential ugi (red) compo- 
nents of the proper motions, together with the line-of-sight 
velocity dispersion a z > (blue). The dispersions are calcu- 
lated along concentric rings from the selected sample of 
2295 stars with proper motions corrected for perspective 
and residual solid-body rotation and 2163 stars with line- 
of-sight velocities corrected for perspective rotation. We 
obtain similar mean velocity dispersion profiles if we only 
use the 718 stars for which both proper motions and line- 
of-sight velocity are measured. We assume the canonical 
distance of 5 kpc to convert the proper motion dispersion 
into units of kms -1 , while the black horizontal lines indi- 
cate the corresponding scale in units of masyr -1 . Below 
the profiles, the crosses represent the corresponding mean 



velocity error per ring, showing that the accuracy of the 
line-of-sight velocity measurements (blue crosses) is about 
four times better than the proper motion measurements 
(green and red crosses, which mostly overlap since the er- 
rors for the two components are similar). 

In § 13.11 we already noticed that since the (smoothed) 
profile of the major-axis proper motion dispersion a x i lies 
on average above that of the minor-axis proper motion 
dispersion ay (Figure and |2J, the velocity distribu- 
tion of u> Cen cannot be fully isotropic. By comparing 
in Figure |S1 the radial (green) with the tangential (red) 
component of the proper motion dispersion, u> Cen seems 
to be radial anisotropic towards the centre, and there is 
an indication of tangential anisotropy in the outer parts. 
Moreover, if lu Cen would be isotropic, the line-of-sight 
velocity dispersion profile (blue) would have to fall on top 
of the proper motion dispersion profiles if scaled with the 
correct distance. A scaling with a distance lower than the 
canonical 5 kpc is needed for the line-of-sight dispersion 
profile to be on average the same as those of both proper 
motion components. 

Hence, it is not surprising that we find a distance as 
low as D = 4.6 ± 0.2 kpc from the ratio of the average 
line-of-sight velocity dispersion and the average proper 
motion dispersion (Appendix Q • This often used simple 
distance estimate is only valid for spherical symmetric ob- 
jects. Whereas the averaged observed flattening for to Cen 
is already as low as q' = 0.879 ± 0.007 (Geyer et al. 1983), 
an inclination of around 48° (§ 14.5(1 . implies an intrinsic 
axisymmctric flattening q < 0.8. 

A model with a constant oblate velocity ellipsoid as in 
Appendix El allows for offsets between the mean velocity 
dispersion profiles. However, the model is not suitable to 
explain the observed variation in anisotropy with radius. 
Therefore, we use in what follows Schwarzschild's method 
to build general axisymmetric anisotropic models. 



5. Schwarzschild's method 

We construct axisymmetric dynamical models using 
Schwarzschild's (1979) orbit superposition method. This 
approach is flexible and efficient and does not require any 
assumptions about the degree of velocity anisotropy. The 
only crucial approximations arc that the object is colli- 
sionless and stationary. While these assumptions are al- 
ways valid for a galaxy, they may not apply to a globular 
cluster. The central relaxation time of u> Cen is a few times 
10 9 years and the half-mass relaxation time a few times 
10 10 years (see also Figure |^ below). The collisionless 
approximation should therefore be fairly accurate. 

The implementation that we use here is an extension 
of the method presented in Verolme et al. (2002). In the 
next subsections, we outline the method and describe the 
extensions. 
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5.1. Mass model 

Schwarzschild's method requires a mass parameterisation, 
which we obtain by using the Multi-Gaussian Expansion 
method (MGE; Monnet, Bacon & Emsellem 1992; 
Emsellem et al. 1994a,b; Cappellari 2002). The MGE- 
mcthod tries to find the collection of two-dimensional 
Gaussians that best reproduces a given surface brightness 
profile or a (set) of images. Typically, of the order of ten 
Gaussians are needed, each with three free parameters: 
the central surface brightness Eo,j, the dispersion along 
the observed major axis o'j and the observed flattening 
q'y Even though Gaussians do not form a complete set of 
functions, in general the surface brightness is well fitted 
(see also Fig. I12|l . Moreover, the MGE-parameterisation 
has the advantage that the deprojection can be performed 
analytically once the viewing angles (in this case the in- 
clination) are given. Also many intrinsic quantities such 
as the potential and accelerations can be calculated by 
means of simple one-dimensional integrals. 

5.2. Gravitational potential 

We deproject the set of best-fitting Gaussians by assum- 
ing that the cluster is axisymmetric and by choosing a 
value of the inclination i. The choice of a distance D to 
the object then allows us to convert angular distances to 
physical units, and luminosities are transformed to masses 
by multiplying with the mass-to-light ratio M/L. 

The latter quantity is often assumed to be independent 
of radius. In the inner regions of most galaxies, where two- 
body relaxation does not play an important role, this often 
is a valid assumption. Generally, globular clusters have 
much shorter relaxation times and may therefore show 
significant M/L-variations. This has been confirmed for 
post core-collapse clusters such as M15 (Dull et al. 1997). 
However, lo Cen has a relatively long relaxation time of 
> 10 9 years, implying that little mass segregation has oc- 
curred and the mass-to-light ratio should be nearly con- 
stant with radius. In our analysis we assume a constant 
M/L, but we also investigate possible M/L-variations. 

The stellar potential is then calculated by applying 
Poisson's equation to the intrinsic density. The contribu- 
tion of a dark object such as a collection of stellar rem- 
nants or a central black hole may be added at this stage. 
On the basis of the relation between the black hole mass 
and the central dispersion (e.g. Tremaine et al. 2002), 
globular clusters might be expected to harbour central 
black holes with intermediate mass of the order 10 3 -10 4 
M Q (e.g. van der Marel 2004). With a central dispersion of 
nearly 20 kms -1 , the expected black hole mass for u> Cen 
would be ~ 10 4 M . The spatial resolution that is required 
to observe the kinematical signature of such a black hole 
is of the order of its radius of influence, which is around 5 
arcsec (at the canonical distance of 5 kpc) . This is approx- 
imately an order of magnitude smaller than the resolution 
of the ground-based observations we use in our analysis, 
so that our measurements are insensitive to such a small 



mass. Hence, we do not include a black hole in our dy- 
namical models of u) Cen. 

5.3. Initial conditions and orbit integration 

After deriving the potential and accelerations, the next 
step is to find the initial conditions for a representative 
orbit library. This orbit library must include all types of 
orbits that the potential can support, to avoid any bias. 
This is done by choosing orbits through their integrals of 
motion, which, in this case, are the orbital energy E, the 
vertical component of the angular momentum L z and the 
effective third integral I3. 

For each energy E, there is one circular orbit in the 
equatorial plane, with radius R c that follows from E = 
$ + \R c d<&/dR c for z = 0, and with §(R,z) the under- 
lying (axisymmetric) potential. We sample the energy by 
choosing the corresponding circular radius R c from a log- 
arithmic grid. The minimum radius of this grid is deter- 
mined by the resolution of the data, while the maximum 
radius is set by the constraint that > 99.9 per cent of the 
model mass should be included in the grid. L z is sampled 
from a linear grid in r\ = £ 2 /L max , with L max the angu- 
lar momentum of the circular orbit. I3 is parameterised 
by the starting angle of the orbit and is sampled linearly 
between zero and the initial position of the so-called thin 
tube orbit (see Figure 3 of Cretton et al. 1999). 

The orbits in the library arc integrated numerically for 
200 times the period of a circular orbit with energy E. In 
order to allow for comparison with the data, the intrin- 
sic density, surface brightness and the three components 
of the projected velocity are stored on grids. During grid 
storage, we exploit the symmetries of the projected veloc- 
ities by folding around the projected axes and store the 
observables only in the positive quadrant (2/ > 0, y' > 0). 
Since the sizes of the polar apertures on which the av- 
erage kinematic data is computed (Figure H3jl are much 
larger than the typical seeing FWHM (1-2 arcsec), we do 
not have to store the orbital properties on an intermedi- 
ate grid and after orbit integration convolve with the point 
spread function (PSF). Instead, the orbital observables are 
stored directly onto the polar apertures. 

5.4. Fitting to the observations 

After orbit integration, the orbital predictions are 
matched to the observational data. We determine the su- 
perposition of orbits whose properties best reproduce the 
observations. If Oij is the contribution of the jih. orbit to 
the ith constraint point, this problem reduces to solving 
for the orbital weights jj in 

No 

Y^'jO., r <- i=l,---,N c , (11) 
3 

where No is the number of orbits in the library, Nq is the 
number of constraints that has to be reproduced and Ci 
is the ith constraint. Since jj determines the mass of each 



G. van de Ven et al.: The dynamical distance and intrinsic structure of u) Cen 



15 



individual orbit in this superposition, it is subject to the 
additional condition 7j > 0. 

Equation (|ll|l can be solved by using linear or 
quadratic programming (e.g. Schwarzschild 1979, 1982, 
1993; Vandervoort 1984; Dejonghe 1989), maximum 
entropy methods (e.g. Richstone & Tremaine 1988; 
Gebhardt et al. 2003) or with a linear least-squares 
solver [e.g. Non-Negative Least-Squares (NNLS), Lawson 
& Hanson 1974], which was used in many of the spherical 
and axisymmetric implementations (e.g. Rix et al. 1997; 
van der Marel et al. 1998; Cretton et al. 1999; Cappellari 
et al. 2002; Verolme et al. 2002; Krajnovic et al. 2005), 
and is also used here. NNLS has the advantage that it is 
guaranteed to find the global best-fitting model and that 
it converges relatively quickly. 

Due to measurement errors, incorrect choices of the 
model parameters and numerical errors, the agreement 
between model and data is never perfect. We therefore 
express the quality of the solution in terms of x 2 j which 
is defined as 



E 



ct-d 

AC, 



(12) 



Here, C* is the model prediction of the constraint C, and 
ACi is the associated error. The value of x 2 f° r a sin- 
gle model is of limited value, since the true number of 
degrees of freedom is generally not known. On the other 
hand, the difference in x 2 between a model and the overall 



minimum value, Ax 



X 



is statistically mean- 



ingful (see Press et al. 1992, § 15.6), and we can assign 
the usual confidence levels to the A\ 2 distribution. The 
probability that a given set of model parameters occurs 
can be measured by calculating A\ 2 for models with dif- 
ferent values of these model parameters. We determine the 
overall best-fitting model by searching through parameter 
space. 

The orbit distribution for the best-fitting model may 
vary rapidly for adjacent orbits, which corresponds to a 
distribution function that is probably not realistic. This 
can be prevented by imposing additional regularisation 
constraints on the orbital weight distribution. This is usu- 
ally done by minimising the nth-order partial derivatives 
of the orbital weights 7j(-E, L z , with respect to the 
integrals of motion E, L z and I3. The degree of smooth- 
ing is determined by the order n and by the maximum 
value A that the derivatives are allowed to have, usually 
referred to as the regularisation error. Since the distribu- 
tion function is well recovered by minimising the second 
order derivatives (n = 2) and smoothening with A = 4 
(e.g. Verolme & de Zeeuw 2002; Krajnovic et al. 2005), we 
adopt these values. 

6. Tests 

Before applying our method to observational data, we test 
it on a theoretical model, the axisymmetric power-law 
model (EZ94). 



6.1. The power-law model 

The potential $ of the power-law model is given by 



$(i?,z) = 



(13) 



in which (R, z) are cylindrical coordinates, $0 is the cen- 
tral potential, R c is the core radius and g$ is the axial 
ratio of the spheroidal equipotentials. The parameter (3 
controls the logarithmic gradient of the rotation curve at 
large radii. 

The mass density that follows from applying Poisson's 
equation to eq. I|13|) can be expanded as a finite sum of 
powers of the cylindrical radius R and the potential 
Such a power-law density implies that the even part of 
the distribution function is a power-law of the two inte- 
grals energy E and angular momentum L z . For the odd 
part of the distribution function, which defines the rota- 
tional properties, a prescription for the stellar streaming 
is needed. We adopt the prescription given in eq. (2.11) of 
EZ94, with a free parameter k controlling the strength of 
the stellar streaming, so that the odd part of the distribu- 
tion function is also a simple power-law of E and L z . Due 
to the simple form of the distribution function, the cal- 
culation of the power-law observables is straightforward. 
Analytical expressions for the surface brightness, the three 
components of the mean velocity and velocity dispersion 
are given in eqs. (3.1)-(3.8) of EZ94. 

6.2. Observables 

We choose the parameters of the power-law model such 
that its observable properties resemble those of u> Cen. 
We use $o=2500 km 2 s -2 , which sets the unit of velocity 
of our models, and a core radius of R c = 2.5 arcmin, which 
sets the unit of length. For the flattening of the potential 
we take g$ = 0.95 and we set (3 = 0.5, so that the even part 
of the distribution function is positive (Fig. 1 of EZ94). 
The requirement that the total distribution function (even 
plus odd part) should be non-negative places an upper 
limit on the (positive) parameter k. This upper limit fc max 
is given by eq. (2.22) of EZ94 7 . Their eq. (2.24) gives the 
value fciso for which the power-law model has a nearly 
isotropic velocity distribution. In our case A; max = 1-38 
and fciso = 1.44. We choose k = 1, i.e., a power-law model 
that has a (tangential) anisotropic velocity distribution. 

Furthermore, we use an inclination of i = 50° , a mass- 
to-light ratio of M/L = 2.5 M Q /L© and a distance of 
D = 5 kpc. At this inclination the projected flattening 
of the potential is q'^ = 0.97. The isocontours of the pro- 
jected surface density are more flattened. Using eq. (2.9) 
of Evans (1994), the central and asymptotic axis ratios of 
the isophotes are respectively q' Q = 0.96 and q'^ — 0.86, 
i.e., bracketing the average observed flattening of w Cen 
of q' = 0.88 (Geyer et al. 1983). 

7 The definition of x has a typographical error and should 
be replaced by x = (1 - 0/2)/|/3|. 
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Fig. 9. Mean velocity and velocity dispersion calculated from a power-law model (first and third column) and from the 
best-fit dynamical Schwarzschild model with D = 4.9 kpc, i — 45° and M/L = 2.5 M q /Lq (second and fourth column). 
The parameters of the power-law model are chosen such that its observables resemble those of u> Cen, including the 
level of noise, which is obtained by randomising the observables according to the uncertainties in the measurements of 
u> Cen (see § 16. 21 and Appendix iBl for details). The average proper motion kinematics in the x'-direction (top row) and 
y'-direction (middle row), and the average mean line-of-sight kinematics (bottom row), calculated in polar apertures 
in the first quadrant, are unfolded to the other three quadrants to facilitate the visualisation. 



Given the above power-law parameters, we calculate 
the three components of the mean velocity V and veloc- 
ity dispersion a on a polar grid of 28 apertures, spanning 
a radial range of 20 arcmin. Because of axisymmetry we 
only need to calculate the observables in one quadrant on 
the plane of the sky, after which we reflect the results to 
the other quadrants. Next, we use the relative precisions 
AV/a ~ 0.11 and Aa/a ~ 0.08 as calculated for u> Cen 
(Appendix [5| , multiplied with the calculated a for the 
power-law model, to attach an error to the power-law ob- 
servables in each aperture. Finally, we perturb the power- 
law observables by adding random Gaussian deviates with 
the corresponding errors as standard deviations. 



Without the latter randomisation, the power-law ob- 
servables are as smooth as those predicted by the dynam- 
ical Schwarzschild models, so that the goodness-of-fit pa- 
rameter x 2 m eq< l|12(l . approaches zero. Such a perfect 
agreement never happens for real data. Including the level 
of noise representative for u> Cen, allows us to use \ 2 to 
not only investigate the recovery of the power-law param- 
eters, but, at the same time, also asses the accuracy with 
which we expect to measure the corresponding parameters 
for uj Cen itself. 



The resulting mean velocity Vobserved and velocity dis- 
persion <7observcd fields for the power-law model are shown 
in respectively the first and third column of Figure[5] They 
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Fig. 10. The (marginalised) goodness-of-fit parameter A\ 2 as a function of distance D, inclination i and mass-to-light 
ratio M/Lv, for different Schwarzschild model fits (indicated by the dots) to an axisymmetric power-law model with 
observables resembling those of oj Cen (see text for details). The x 2 -values are offset such that the overall minimum, 
indicated by the cross, is zero. The contours are drawn at the confidence levels for a A% 2 -distribution with three 
degrees of freedom, with inner three contours corresponding to the 68.3%, 95.4% and 99.7% (thick contour) confidence 
levels. Subsequent contours correspond to a factor of two increase in A% 2 . The input parameters D = 5.0 kpc, i = 50° 
and M/L = 2.5 M /L , indicated by the open square, are recovered within the 68.3% confidence levels. 



are unfolded to the other three quadrants to facilitate the 
visualisation. 

6.3. Schwarzschild models 

We construct axisymmetric Schwarzschild models based 
on the power-law potential l|13|l . We calculate a library 
of 2058 orbits by sampling 21 energies E, 14 angular mo- 
menta L z and 7 third integrals I3 . In this way, the number 
and variety of the library of orbits is large enough to be 
representative for a broad range of stellar systems, and 
the set of eqs. (f 1 1 |l is still solvable on a machine with 512 
Mb memory (including regularisation constraints). 

The resulting three-integral Schwarzschild models in- 
clude the special case of dependence on only E and L z 
like for the power-law models. Schwarzschild's method re- 
quires that the orbits in the library are sampled over a 
range that includes most of the total mass, whereas all 
power-law models have infinite mass. To solve this prob- 
lem at least partially, we ensure that there are enough 
orbits to constrain the observables at all apertures. We 
distribute the orbits logarithmically over a radial range 
from 0.01 to 100 arcmin (five times the outermost aper- 
ture radius) and fit the intrinsic density out to a radius of 
10 5 arcmin. The orbital velocities are binned in histograms 
with 150 bins, at a velocity resolution of 2 kms -1 . 

To test whether and with what precision we can re- 
cover the input distance of D = 5 kpc, the inclination of 
i = 50° and the mass-to-light ratio M/L = 2.5 Mq/L , 
we calculate models for values of D between 3.5 and 6.5 
kpc, i between 35° (the asymptotic isophotal axis ratio 
q'^ = 0.86 implies that i > 30°) and 70°, and M/L 
between 1.5 and 3.5 Mq/L q . Additionally, to test how 
strongly the best-fitting parameters depend on the un- 



derlying mass model, we also vary the flattening of the 
power-law potential (7$ between 0.90 and 1.00. 

We then fit each of the dynamical models simultane- 
ously to the calculated observables of the power-law model 
(with q<£ = 0.95). Comparing these calculated observ- 
ables with those predicted by the Schwarzschild models, 
results for each fitted Schwarzschild model in a goodness- 
of-fit parameter % 2 . We use this value to find the best-fit 
Schwarzschild model and to determine the accuracy of the 
corresponding best-fit parameters. 

Calculating the observables for all orbits in the library 
requires about an hour on a 1 GHz machine with 512 MB 
memory and the NNLS-fit takes about half an hour. No 
distinct models need to be calculated for different values 
of M/L, as a simple velocity scaling prior to the NNLS-fit 
is sufficient. Making use of (a cluster of) about 30 com- 
puters, the calculations for the full four-parameter grid of 
Schwarzschild models takes a few days. 

6.4. Distance, inclination and mass-to-light ratio 

The Schwarzschild model that best fits the calculated 
power-law observables is the one with the (overall) low- 
est x 2 -value. After subtraction of this minimum value, 
we obtain A\ 2 as function of the three parameters D, 
i and M/L (with g$ = 0.95 fixed). To visualise this three- 
dimensional function, we calculate for a pair of parame- 
ters, say D and i, the minimum in A% 2 as function of 
the remaining parameter, M/L in this case. The contour 
plot of the resulting marginalised A% 2 is shown in the left 
panel of Figure ^] The dots indicate the values at which 
Schwarzschild models have been constructed and fitted 
to the power-law observables. The contours are drawn at 
the confidence levels for a Ax 2 -distribution with three de- 
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Fig. 11. The (marginalised) goodness-of-fit parameter A% 2 as a function of distance D, inclination i and mass-to-light 
ratio M/L against the flattening q§ of the underlying potential, for different Schwarzschild model fits (indicated by 
the dots) to the observables of an axisymmetric power-law model resembling those of ui Cen. The contours are as 
in Figure ITU1 but for a A% 2 -distribution with four degrees of freedom. The cross indicates the overall best-fit model 
(A% 2 = 0). The input parameters of the power-law model, q$> = 0.95, D = 5.0 kpc, i = 50° and M/L = 2.5 M Q /L , 
are indicated by the open square. The input parameters are recovered within the 68.3% confidence levels, even for mass 
models that assume a (slightly) incorrect value for the flattening. However, spherical models (g$ = 1.0) are strongly 
ruled out. 



grees of freedom, with inner three contours corresponding 
to the 68.3%, 95.4% and 99.7% (thick contour) confidence 
levels. Subsequent contours correspond to a factor of two 
increase in A% 2 . The minimum (A% 2 = 0) is indicated by 
the cross. Similarly, we show in the middle and left panel 
the contour plots of A\ 2 marginalised for respectively the 
pair D and M/L and the pair i and M/L. 

The input parameters D — 5.0 kpc, i = 50° and 
M/L = 2.5 M Q /L Q , indicated by the open square, are well 
recovered. The mean velocity Vmodei and velocity disper- 
sion <7 mo dei predicted by the best-fit Schwarzschild model 
are shown in the second and fourth column of Figure 
The corresponding power-law observables are well repro- 
duced within the error bars. 

Since the parameters of the power-law model are cho- 
sen such that its observables and corresponding errors re- 
semble those of u> Cen, the contours in Figure ^| pro- 
vide an estimate of the precision with which we expect 
to measure the best-fitting parameters for u) Cen. At the 
68.3%-level (99.7%-level) the distance D, inclination i and 
mass-to-light ratio M/L are retrieved with an accuracy of 
respectively 6 (11), 9 (18), 13 (28) per cent. Due the ad- 
ditional complication of infinite mass in the case of the 
power-law models, these estimates most likely are upper 
limits to the precision we expect to achieve for to Cen. This 
holds especially for the inclination and the mass-to-light 
ratio as they are sensitive to how well the mass model 
is fitted. The distance is mainly constrained by the kine- 
matics, so that the corresponding accuracy is probably an 
accurate estimate of the precision with which we expect 
to measure the distance to u> Cen. 



6.5. Flattening 

The above investigation of the recovery of the global pa- 
rameters D, i and M/L is for a known mass model, given 
by the power-law potential (|13|l . In general, we obtain the 
mass model from a MGE-parameterisation of the observed 
surface brightness (§ IS.lfl . There is no guarantee that the 
resulting MGE model provides an accurate description of 
the true mass distribution. Therefore, we tested the ef- 
fect of an incorrect mass model on the best-fit parameters 
by varying the flattening g$ of the power-law potential 
while keeping the calculated observables (for the power- 
law model with q$ = 0.95) fixed. 

Since we use these same values for the other power-law 
parameters ($o=2500 km 2 s~ 2 , R c = 2.5 arcmin, (3 = 0.5 
and k = 1), we have to be careful that by varying q$ the 
resulting model is still physical, i.e., that the underlying 
distribution function is non-negative. For these parame- 
ters and <?$ between 0.9 and 1.0 this is the case (EZ94). 

As before, for all Schwarzschild models we calculate 
A% 2 , which is now a function of the four parameters D, 
i, M/L and <?$. In the three panels of Figure ITTI we show 
A% 2 marginalised for respectively D, i and M/L against 
q$ . The symbols and contours are as in Figure 1101 but 
now for a Ax 2 -distribution with four degrees of freedom. 
The input parameters of the power-law model (indicated 
by an open square) are g$ = 0.95, D = 5.0 kpc, i = 50° 
and M/L = 2.5 M /L Q . 

The distance D is well constrained around the correct 
input value, even at values for the flattening of the po- 
tential g$ that arc different from the true value of 0.95. 
This implies that the best-fitting distance is accurate even 
for mass models that assume a (slightly) incorrect value 
for the flattening. Whereas a potential with a flattening 
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Table 2. The parameters of the 8 Gaussians from the 
MGE-fit to the V-band surface brightness profile of u> Cen 
as derived by Meylan (1987). The second column gives the 
central surface brightness (in L Q pc~ 2 ) of each Gaussian 
component, the third column the dispersion (in arcmin) 
along the major axis and the fourth column the observed 
flattening. 



3 




a 


q 




(L© PC" 2 ) 


(arcmin) 




1 


2284.7077 


0.15311 


1.000000 


2 


3583.7901 


1.47715 


0.934102 


3 


3143.2029 


2.52542 


0.876713 


4 


1670.8477 


3.69059 


0.848062 


5 


840.86244 


5.21905 


0.849760 


6 


262.69433 


7.53405 


0.835647 


7 


46.995893 


11.0685 


0.866259 


8 


3.3583961 


17.5470 


0.926328 



as low as 0.90 still (just) falls within the contour at the 
99.7%-level, we conclude, as in § 14.61 that spherical mod- 
els (g$ = 1) are strongly ruled out. The middle and right 
panel of Figure ^2 show that the results for respectively 
the mass-to-light M/L and inclination i are similar, al- 
though, as before, they are less well constrained due to 
the infinite mass of the power-law models. 

7. Dynamical models for to Cen 

We use our method as described in § \5\ to construct dy- 
namical models for lu Cen. We obtain the mass model from 
a MGE-parametrisation of the observed surface bright- 
ness. We compute the mean velocity and velocity disper- 
sion of both proper motion components and along the line- 
of-sight in polar apertures on the plane of the sky. For 
a range of distances, inclinations and (constant) mass- 
to-light ratios, we then simultaneously fit axisymmetric 
Schwarzschild models to these observations. Additionally, 
we also allow for radial variation in the mass-to- light ratio. 

7.1. MGE mass model 

An MGE-fit is best obtained from a two-dimensional im- 
age, which gives direct information about the flattening 
and any radial variations in the two-dimensional structure 
of the object. Unfortunately, no such image is available to 
us, and the only published surface brightness observations 
of ui Cen consist of radial surface brightness profiles, and 
an ellipticity profile by Geyer et al. (1983). We therefore 
perform a one-dimensional MGE-fit to the radial surface 
brightness profile, and after that use the ellipticity profile 
to include flattening in the mass model. 

We use the U-band surface brightness data from 
Meylan (1987), who combined various published measure- 
ments (Gascoigne & Burr 1956; Da Costa 1979; King et 
al. 1968). Their data consists of individual measurements 
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Fig. 12. The Multi-Gaussian Expansion (MGE) of the V- 
band surface brightness profile of u> Cen. The filled circles 
represent the measurements by Meylan (1987), the dotted 
curves correspond to the eight Gaussians in the expansion 
and the solid curve represents their sum. The top panel 
shows the surface brightness £ as a function of projected 
radius R' (in arcmin). Kalnajs (1999) has shown that the 
quantity R'Y, in the bottom panel is a good diagnostic of 
the mass that is enclosed at each radius. 



along concentric rings, while the MGE-algorithm devel- 
oped by Cappellari (2002) requires a regular (logarith- 
mic) spacing of the surface brightness measurements. We 
therefore first describe the profile in terms of a fourth- 
order polynomial and then fit a set of one-dimensional 
Gaussians to this polynomial. Eight Gaussians with dif- 
ferent central surface brightness Eqvj and dispersion a'j 
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Table 3. The mean velocity and velocity dispersion calculated in polar apertures on the plane of sky from the proper 
motion observations. Per row the information per aperture is given. The first column labels the aperture and the 
second column gives the number of stars n* that fall in the aperture. Columns 3-6 list the polar coordinates r (in 
arcmin) and the angle 9 (in degrees) of the centroid of the aperture and the corresponding widths Ar (in arcmin) and 
AO (in degrees). The remaining columns present the average proper motion kinematics in units of masyr -1 . The mean 
velocity V with error AV and velocity dispersion a with error Act are given in columns 7-10 for the proper motion 
component in the x'-direction and in columns 11-14 for the proper motion component in the y'-direction. 



1 

2 
3 
4 
5 
6 
7 
8 
9 
10 
11 
12 
13 
14 
15 
16 
17 
18 
19 
20 
21 
22 
23 
24 
25 
26 
27 
28 



80 
90 
67 
74 
85 
77 
76 
82 
105 
88 
70 
72 
65 
95 
88 
91 
73 
72 
61 
88 
95 
61 
85 
68 
58 
74 
79 
92 



ro 
1.14 
3.04 
3.04 
3.04 
4.59 
4.59 
4.59 
4.59 
6.31 
6.31 
6.31 
6.31 
6.31 
8.49 
8.49 
8.49 
8.49 
8.49 
8.49 
11.54 
11.54 
11.54 
11.54 
11.54 
16.64 
16.64 
16.64 
16.64 



45.0 
15.0 
45.0 
75.0 
11.2 
33.7 
56.2 
78.7 

9.0 
27.0 
45.0 
63.0 
81.0 

7.5 
22.5 
37.5 
52.5 
67.5 
82.5 

9.0 
27.0 
45.0 
63.0 
81.0 
11.2 
33.7 
56.2 
78.7 



Ar 
2.28 
1.53 
1.53 
1.53 
1.57 
1.57 
1.57 
1.57 
1.86 
1.86 
1.86 
1.86 
1.86 
2.52 
2.52 
2.52 
2.52 
2.52 
2.52 
3.56 
3.56 
3.56 
3.56 
3.56 
6.64 
6.64 
6.64 
6.64 



A6> 
90.0 
30.0 
30.0 
30.0 
22.5 
22.5 
22.5 
22.5 
18.0 
18.0 
18.0 
18.0 
18.0 
15.0 
15.0 
15.0 
15.0 
15.0 
15.0 
18.0 
18.0 
18.0 
18.0 
18.0 
22.5 
22.5 
22.5 
22.5 



-0.15 
-0.16 

0.03 
-0.15 
-0.27 
-0.08 
-0.20 
-0.19 

0.00 
-0.13 
-0.28 
-0.25 
-0.25 
-0.04 
-0.09 
-0.15 
-0.31 
-0.35 
-0.40 

0.02 
-0.17 
-0.24 
-0.41 
-0.36 
-0.02 
-0.14 
-0.17 
-0.21 



0.09 
0.07 
0.12 
0.08 
0.06 
0.07 
0.07 
0.05 
0.06 
0.07 
0.07 
0.05 
0.07 
0.05 
0.05 
0.05 
0.06 
0.05 
0.07 
0.05 
0.04 
0.05 
0.05 
0.05 
0.06 
0.06 
0.05 
0.05 



0.80 
0.66 
0.90 
0.64 
0.57 
0.63 
0.55 
0.55 
0.60 
0.61 
0.58 
0.45 
0.58 
0.56 
0.46 
0.49 
0.51 
0.44 
0.58 
0.44 
0.42 
0.44 
0.44 
0.43 
0.40 
0.48 
0.46 
0.43 



0.07 
0.04 
0.07 
0.07 
0.03 
0.05 
0.05 
0.04 
0.04 
0.04 
0.07 
0.04 
0.05 
0.04 
0.04 
0.04 
0.06 
0.04 
0.05 
0.04 
0.04 
0.04 
0.03 
0.03 
0.04 
0.05 
0.03 
0.03 



Yvl 
-0.01 

0.23 

0.06 

-0.08 
0.19 
0.13 
0.13 
0.07 
0.26 
0.13 
0.23 

-0.01 
0.05 
0.22 
0.10 
0.14 
0.19 
0.14 

-0.03 
0.20 
0.17 
0.18 
0.05 
0.05 
0.19 

-0.01 
0.04 

-0.05 



AVy, 

0.09 
0.07 
0.08 
0.09 
0.06 
0.06 
0.08 
0.07 
0.05 
0.05 
0.06 
0.06 
0.06 
0.04 
0.07 
0.04 
0.05 
0.06 
0.07 
0.05 
0.05 
0.05 
0.04 
0.05 
0.06 
0.06 
0.04 
0.04 



0.70 
0.64 
0.62 
0.71 
0.57 
0.57 
0.69 
0.66 
0.50 
0.48 
0.50 
0.53 
0.45 
0.38 
0.53 
0.41 
0.44 
0.54 
0.48 
0.46 
0.49 
0.41 
0.43 
0.46 
0.41 
0.45 
0.41 
0.35 



Acy 
0.05 
0.05 
0.05 
0.06 
0.05 
0.08 
0.06 
0.06 
0.04 
0.05 
0.06 
0.05 
0.03 
0.02 
0.07 
0.03 
0.03 
0.05 
0.04 
0.04 
0.04 
0.03 
0.03 
0.03 
0.05 
0.04 
0.04 
0.03 



are required by the MGE-fit (second and third column of 
Table [5J. Figure IT21 shows that this MGE- model provides 
an excellent fit, not only to the surface brightness S, but 
also to R'T, (cf. Kalnajs 1999). 

The MGE-parameterisation is converted into a two- 
dimensional luminosity distribution by assigning an ob- 
served flattening q'j to each Gaussian in the superposi- 
tion. We take into account that the observed flattening 
of to Cen varies as a function of radius (cf. Geyer et 
al. 1983). This is done by assuming that the flattening 
of the jth Gaussian q'j is equal to the observed flatten- 
ing at a projected radius R' = a'j. This is justified by 
the fact that a given Gaussian contributes most at radii 
close to its dispersion cr'-. Although small deviations from 
the true two-dimensional light distribution in to Cen may 
still occur, we showed in § 16.51 that this approximation 
does not significantly influence the derived intrinsic pa- 
rameters for to Cen. Moreover, a two-dimensional MGE-fit 
to the combination of the surface brightness profile from 
Meylan (1987) and the ellipticity profile from Geyer et 



al. (1983), yields nearly equivalent MGE parameters as 
those in Table |2 although the fit to the observed surface 
brightness profile is less good. 

To conserve the total luminosity, we increase the cen- 
tral surface brightness of each Gaussian by l/q'j- Taking 
into account a reddening of E(B — V) = 0.11 for to Cen 
(Lub 2000), the total I^-band luminosity of our mass 
model, at the canonical distance of 5.0 ± 0.2 kpc, is 
Ly — 1.0 ± 0.1 x 10 6 L©. This compares well with other 
estimates of the total luminosity of to Cen of 0.8 x 10 6 
L Q (Carraro & Lia 2000), 1.1 x 10 6 L Q (Seitzer 1983) and 
1.3 x 10 6 L (Meylan 1987). The most flattened Gaussian 
in the superposition (J = 7) places a mathematical lower 
limit on the inclination of 33°. This is safely below the 
constraint 41-57 degrees found in § 14.51 

7.2. Mean velocity and velocity dispersion 

We construct a polar aperture grid for the proper motions 
and line-of-sight velocities, as shown in Figure ED The 
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dots in the top panel represent the positions, folded to 
the first quadrant, of the 2295 selected stars with ground- 
based proper motions. The overlayed polar grid, extending 
to about 20 arcmin, consists of 28 apertures. Per aperture, 
the number of stars is indicated, adding up to a total of 
2223 stars. Similarly, the bottom panel shows the 2163 
selected stars with line-of-sight velocities. The different 
number of stars and spatial distribution results in a polar 
grid of 27 apertures, which includes in total 2121 stars. 

For each aperture, we use the maximum likelihood 
method (Appendix [SJ to compute the mean velocity V 
and velocity dispersion a for both proper motion compo- 
nents on along the linc-of-sight. We calculate correspond- 
ing errors by means of the Monte Carlo bootstrap method. 

Each aperture contains around 50 to 100 stars. In 
Appendix iBl we find that this is a good compromise be- 
tween precision in the observables and spatial resolution. 
Including more stars per aperture by increasing its size 
decreases the uncertainties in the observables (and hence 
makes the resulting kinematic fields smoother). At the 
same time, since the apertures should not overlap to as- 
sure uncorrelated observables, this means less apertures 
in the polar grid and hence a loss in spatial resolution. 

The properties of the apertures and corresponding 
mean kinematics are given in Table for the proper mo- 
tions and in Table 0] for the line-of-sight velocities. The 
mean velocity Vobscrved and velocity dispersion <7 bservcd 
fields are shown in the first and third column of Figure ITU 
respectively. Although the average kinematics are only cal- 
culated in the first quadrant, we use the assumed axisym- 
metric geometry to unfold them to the other three quad- 
rants to facilitate the visualisation. 

7.3. Constructing dynamical models 

First, we calculate models for a range of values in distance 
D, inclination i and constant y-band mass-to-light ratio 
M/Ly- Next, fixing D and i at their measured best-fit 
values, we also calculate a large set of models in which we 
allow M/Ly to vary with radius. 

We sample the orbits on a grid of 21 x 14 x 7 values 
in (E,L z ,l3) on a radial range from 0.01 to 63 arcmin. 
This grid extends beyond the tidal radius of 45 arcmin 
Trager et al. 1995), so that all mass is included. No PSF- 
convolution is used and the observables are stored directly 
onto the apertures. 

We (linearly) sample D between 3.5 and 6.5 kpc in 
steps of 0.5 kpc, and additionally we refine the grid be- 
tween 4.0 and 5.5 kpc to steps of 0.1 kpc. We vary i be- 
tween 35 (close to the lower limit of 33 degrees imposed 
by the flattening, see § 17. 1[) and 90 degrees in steps of five 
degrees, and we refine between 40 and 50 degrees to steps 
of one degree. We choose the constant M/Ly values be- 
tween 2.0 and 4.0 M Q /L with steps 0.5 M0/L0, and we 
refine between 2.0 and 3.0 M Q /L to steps of 0.1 Mq/Lq. 

To investigate possible variation in M/Ly with radius, 
we make use of the eight Gaussian components of the 



Table 4. The mean velocity and velocity dispersion cal- 
culated in polar apertures on the plane of sky from the 
line-of-sight velocity observations. Columns 1-6 are as in 
Table El and the remaining columns present the average 
line-of-sight kinematics in kms . 
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MGE mass model (§ 17. 1|) . In case of constant M/Ly, we 
obtain the intrinsic density by multiplying all the (depro- 
jected) components with the same constant M/Ly value. 
To construct a mass model with a radial M/Ly profile, 
we multiply each component with its own M/Ly value, as 
in this way the calculation of the potential is still efficient. 

However, to reduce the number of free parameters (to 
make a search through parameter space feasible) and to 
enforce a continuous profile, we only vary the M/Ly val- 
ues for the first, second, fourth and sixth component. For 
the third and fifth component, we interpolate between 
the M/L values of the neighbouring components. To the 
outer two components we assign the same M/Ly value 
as the sixth component, because their individual M/Ly 
values are not well constrained due to the small number 
of kinematic measurements at these radii. With the dis- 
tance and inclination fixed at their best-fit values from 
the case of constant mass-to-light ratio, we are left with a 
four-dimensional space to search through, requiring again 
a few days on (a cluster of) about 30 computers. 

All dynamical models are fitted simultaneously to the 
two-dimensional light distribution of u> Cen (S I7.1|I . and to 
the mean velocity and velocity dispersion of both proper 
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Fig. 14. Mean velocity and velocity dispersion calculated from the observations of uj Cen (first and third column) 
and from the best-fit dynamical model with D = 4.8 kpc, i = 50° and MjLy = 2.5 Mq/Lq (second and fourth 
column). The mean proper motion kinematics in the aZ-direction (top row) and y'-direction (middle row), and the 
mean line-of-sight kinematics (bottom row), calculated in polar apertures in the first quadrant, are unfolded to the 
other three quadrants to facilitate the visualisation. 



motions components and along the line-of-sight, calcu- 
lated in polar apertures on the plane of the sky (first and 
third column of Figure [Tty . Comparing the predicted val- 
ues with the observations, results for each fitted model in 
a goodness-of-fit parameter x 2 . We use this value to find 
the best-fit model and to determine the accuracy of the 
corresponding best-fit parameters. 

8. Best-fit parameters 

In Figure El we show A% 2 as a (marginalised) function 
of the distance D, inclination i and constant mass-to- 
light ratio MjLy. The dots represent the values at which 
dynamical models have been constructed and fitted to 
the two-dimensional (photometric and kinematic) obser- 
vations of ui Cen. The cross indicates the over-all best-fit 



model. The contours show that all three parameters are 
tightly constrained, with at the 68.3%-level (99.7%-level): 
D = 4.8 ± 0.3 (±0.5) kpc, % = 50 ± 3 (±5) degrees and 
MjLy = 2.5 ±0.1 (±0.2) M Q /L Q . As an illustration that 
our best-fit model indeed reproduces the observations, the 
mean velocity and velocity dispersion in polar apertures 
on the plane of the sky as they follow from this model 
are shown in respectively the second and fourth column 
of Figure ^] The model fits the observations within the 
uncertainties given in Table and [JJ 

After the discussion on the set of models where we 
allow the mass-to-light ratio MjLy to vary with radius, 
we compare our best-fit values for the (constant) mass- 
to-light ratio, inclination and distance with results from 
previous studies. 
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Fig. 15. The (marginalised) goodness-of-fit parameter A\ 2 as a function of distance D, inclination i and mass-to-light 
ratio M/Lv, for different dynamical model fits (indicated by the dots) to the kinematics of u> Cen. The contours are 
as in Figure ITU1 The best-fit dynamical model is at D = 4.8 kpc, i = 50° and M/Lv = 2.5 M Q /L , indicated by the 
cross The dashed curve shows the Dtani = 5.6 kpc constraint from the mean velocities f§ 14. 5|) . 



8.1. Mass-to-light ratio variation 

Figure 1161 summarises the results from fitting models in 
which we allowed the mass-to-light ratio M/Lv to vary 
with radius in the way described in § 17.31 The filled cir- 
cles represent the eight Gaussian components, with the 
best-fit M/Lv value of each component plotted against 
their dispersions along the major axis (see column three 
of TableEJ- The error bars represent the 68.3% confidence 
level. 

The uncertainty on the innermost point around 0.15 
arcmin is relatively large since at that small radius there 
are only a few observations (see Figure to constrain 
the M/Ly value. Nevertheless, the resulting M/Lv profile 
only shows a small variation, which is not significantly 
different from the best-fit constant M/Lv of 2.5 M Q . 

In the above experiment, we fixed the distance and in- 
clination at the best-fit values of D = 4.8 kpc and i = 50 
degrees from the case of constant M/Lv- Although an im- 
portant constraint is that all eight Gaussian components 
have to be at the same distance, its precise value, as well 
as that of the inclination, is not crucial. We tested that 
a reasonable variation in these fixed values (within the 
99.7% confidence level in Figure ITS)) does not significantly 
change the best-fit M/Lv profile. We conclude that a con- 
stant mass-to-light ratio for to Cen is a valid assumption. 

8.2. Mass-to-light ratio 

Our best-fit mass-to-light ratio of M/Lv — 2.5 ± 0.1 
Mq/Lq lies in between the estimates by Seitzer (1983) 
of 2.3 Mq/Lq and by Meylan (1987) of 2.9 Mq/L . 
Meylan et al. (1995) derived a value of 4.1 Mq/Lq, based 
on a spherical, radial anisotropic King-Michie dynamical 
model, while we find that u> Cen is flattened and out- 
wards tangentially anisotropic (see § 19.2(1 . Moreover, their 
adopted central value of the line-of-sight velocity disper- 



sion is significantly higher than ours, even if we use the 
same data-set by M97. 

Meylan et al. (1995) estimated the total mass of to Cen 
to be 5.1 x 10 6 M Q , which is also significantly higher than 
what we derive. After multiplication with the total lumi- 
nosity of our mass model of L = 1.0 x 10 6 L Q (at the 
best-fit distance of D = 4.8 ± 0.3 kpc), we find a total 
mass of M = (2.5 ±0.3) x 10 6 M Q . This is consistent with 
the value by Mandushev et al. (1991) of 2.4 x 10 6 M Q and 
Seitzer (1983) of 2.8 x 10 6 M©. The estimate by Meylan 
(1987) of 3.9 x 10 6 M Q is higher, but again based on a 
spherical King-Michie model. 

8.3. Inclination 

The dashed curve in the left panel of Figure El shows 
the Dt'Aiii = 5.6 kpc constraint from the mean velocities 
derived in § 14.51 This constraint can be used to eliminate 
either the distance or the inclination and hence reduce the 
parameter space. Although we do not use this constraint 
in the dynamical models, it is clear that the above best-fit 
D and i yield Dtani = 5.6 ± 0.2 kpc, which is consistent 
with the value derived from the mean velocities. 

The best-fit inclination of i ~ 50 ± 3 degrees falls 
within the range of 30-60 degrees that was derived in 
Paper I from the amplitude of the proper motions, but 
is slightly higher than the estimate by van Leeuwen & Le 
Poole (2002) between 40 and 60 degrees. However, as dis- 
cussed in § 14.51 they used models of modest complexity 
and freedom which require strong assumptions, whereas 
our method is much more general and robust. 

Our best-fit inclination implies that u> Cen is intrinsi- 
cally even more non-spherical than the average observed 
flattening of q' = 0.879 ±0.007 (Geyer et al. 1983) already 
indicates. Using the relation q 2 sin 2 i = q' 2 — cos 2 i for ax- 
isymmetric objects, we find an average intrinsic axial ratio 
q = 0.78 ±0.03. 
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Fig. 13. The polar aperture grid for the proper motions 
(top panel) and for the line-of-sight velocities (bottom 
panel). The dots represent the individual stars, with posi- 
tions folded to the first quadrant, while the solid lines in- 
dicate the locations of the apertures. The number of stars 
included are indicated in each aperture. An enlargement 
of the inner part of the line-of-sight polar grid is shown in 
the top-right corner of the bottom panel. 

8.4. Distance 

Adopting a reddening of E(B — V) = 0.11 for uo Cen (Lub 
2000), the best-fit dynamical distance corresponds to a 
distance modulus of (to - M) v = 13.75 ± 0.13 (±0.22 at 
the 99.7%- level). This is consistent with the (canonical) 
distance modulus of (to — M)y = 13.84 by photomet- 
ric methods, as given in the globular cluster catalogue of 
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Fig. 16. Variation in mass-to-light ratio M/Ly with pro- 
jected radius R' . The filled circles represent the eight 
Gaussian components of the MGE mass model, with the 
best-fit M/Ly value of each component plotted against 
its dispersion along the major axis. With the distance and 
inclination fixed at D = 4.8 kpc and i = 50 degrees, we 
allowed variation in the M/Ly values for the four inner 
points with error bar, while the two outer points were 
shifted vertically similar to the fourth point, and the re- 
maining two points were interpolated between the two 
neighbouring points. Each of the models was simultane- 
ously fitted to the photometric and kinematic observations 
of uj Cen. The error bars represent the 68.3% confidence 
level for the corresponding Ax 2 -distribution with four de- 
grees of freedom. The variation in the resulting Mj Ly pro- 
file is small with no significant deviation from the best-fit 
constant M/Ly of 2.5 M Q (horizontal dashed line). 

Harris et al. (1996), together with the uncertainty estimate 
of about 0.1 magnitude by Benedict et al. (2002), using the 
absolute magnitude of RR Lyrae stars. Using the infrared 
colour versus surface brightness relation for the eclipsing 
binary OGLEC 17, Thompson et al. (2001) find a larger 
distance modulus of (to — M)y = 14.05 ± 0.11. However, 
their distance modulus estimates based on the measured 
bolometric luminosity of the binary components, are on 
average lower, ranging from 13.66 to 14.06. 

Although our dynamical distance estimate is consis- 
tent with that by other methods, it is at the lower end. A 
lower value for the distance is expected if the proper mo- 
tion dispersion is over-estimated and/or the line-of-sight 
velocity dispersion under-estimated (see also Appendix IO 
eq. IC.lll . As we saw in § EI both are likely in the case of 
uj Cen if the kinematic data is not properly selected. The 
correction in § 0] for perspective rotation and especially 
for the residual solid-body rotation is crucial for the con- 
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struction of a realistic dynamical model and also to find a 
reliable distance estimate. 

An impression of the effect of the selection and correc- 
tion of the kinematic data on the distance estimate follows 
from the range of dynamical models we constructed for 
u> Cen. Before any selection and correction, the kinematics 
of the cluster stars give rise to a best-fit dynamical model 
at a distance as low as ^3.5 kpc. After removing from 
the proper motion data-set the stars disturbed by their 
neighbours, i.e., only selecting class stars, the best-fit 
distance becomes ^4.0 kpc. The correction for perspec- 
tive and solid-body rotation increase the best-fit distance 
to ^4.5 kpc. Finally, after the additional selection on ve- 
locity errors, we find our best-fit dynamical distance of 
4.8 ±0.3 kpc. 

An even tighter selection does not significantly change 
the best-fit dynamical model and corresponding distance. 
The same is true if we use a different polar grid, with 
fewer or more stars per aperture, and if we restrict to 
only fitting the average kinematics in the inner or outer 
parts. Still, e.g. remaining interlopers in the proper motion 
data-set can cause a (small) under-estimation of the dis- 
tance. Moreover, Platais et al. (2003) argue that possibly 
a (non-physical) residual proper motion colour/magnitude 
dependence in the data-set of Paper I causes the system- 
atic offset between the proper motions of the metal-rich 
RGB-a stars and those of the dominant HB and metal- 
poor RGB stars, noticed by Ferraro, Bellazzini & Pancino 
(2002). Since we do not correct for this possible system- 
atic offset, the proper motion dispersion might be over- 
estimated and hence our distance estimate can be sys- 
tematically too low. However, the effect is expected to be 
small since the number of RGB-a stars in the data-set is 
small. A deeper proper motion catalogue, like that of King 
& Anderson 2002) obtained with the HST, is needed to 
better quantify (non-physical and physical) differences in 
the proper motions among the multiple stellar populations 
observed in w Cen. 

Although the distance and inclination are tightly 
linked through the mean velocities (§ 14.5(1 , a small under- 
estimation of the distance only results in a slight over- 
estimation of the inclination (see also the solid curve in 
the right panel of Figure 0) . Similarly, the mass-to-light 
ratio is nearly insensitive to small changes in the distance. 

9. Intrinsic structure 

We use the intrinsic velocity moments of our best-fit dy- 
namical model to investigate the importance of rotation 
and the degree of anisotropy in u> Cen. Additionally, the 
distribution of the orbital weights allows us to study the 
phase-space distribution function of ui Cen. 

9.1. Rotation 

We calculate the intrinsic velocity moments of our best-fit 
model by combining the appropriate moments of the or- 
bits that receive weight in the superposition. We consider 
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Fig. 17. The colours represent the mean azimuthal rota- 
tion (v^) in the meridional plane as a function of equa- 
torial plane radius R and height z, and normalised by 
crms (excluding the axes to avoid numerical problems). 
The black curves are contours of constant mass density in 
steps of one magnitude, from the mass model (solid) and 
from the best-fit model (dashed), showing that the mass 
is well fitted. 



the first and second order velocity moments, for which 
(vr) = (ve) = (vRVj,) — (vgVff,) — because of axisym- 
mctry. We define the radial, angular and azimuthal ve- 
locity dispersion respectively as o\ — (v R ), <Jg = (v 2 ), 
°4> ~ ( v 4>) — ( v <t>) • ^ ne on ly non- vanishing cross- term is 
a R8 — ( V R V 4>)- The total root-mean-square velocity dis- 
persion ctrms is given by cr RMS = (a% + <r 2 + cr£)/3. 

A common way to establish the importance of rota- 
tion in elliptical galaxies and bulges of disk galaxies, is 
to determine their position in the (V/cr, e)-diagram (e.g. 
Davies et al. 1983). The observational quantities that are 
used for V, a and e are respectively the maximum (line-of- 
sight) velocity along the major axis, the average velocity 
dispersion within half the effective radius and the ellip- 
ticity at the effective radius. We obtain for ui Cen the 
observational quantities V ~ 8 kms -1 (at a radius of ~ 8 
arcmin), a ~ 16 kms -1 and e ~ 0.15 (Geyer et al. 1983). 
These values result in (V/cr, e) ~ (0.5, 0.15), placing w Cen 
just above the curve for isotropic oblate rotators. 

On the other hand, the intrinsic velocity moments from 
our best-fit dynamical model for uj Cen, allow us to investi- 
gate intrinsically the importance of rotation. The colours 
in Figure 1171 show the ratio of the mean (azimuthal) ro- 
tation (vq) over the total root-mean-square velocity dis- 
persion crms j a s function of the position in the merid- 
ional plane. Near the equatorial plane and between radii 
of about 5 to 15 arcmin, this ratio is > 0.5. The maxi- 
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mum of ~ 0.7 around 8 arcmin coincides with the peak 
in the mean line-of-sight velocity field. Within this region 
in the meridional plane rotational support is important. 
However, more inwards and further outwards this ratio 
rapidly drops below 0.5 and u> Cen is at least partly pres- 
sure supported. We conclude that rotation is important in 
u> Cen, but it is not a simple isotropic oblate rotator. 

9.2. Anisotropy 

For the velocity distribution in lu Cen to be isotropic all 
three velocity dispersion components an, og and cr^ have 
to be equal and the cross-term a Re has to vanish. Figure lTFI 
shows that this is not the case. 

In the upper panels, we show the degree of anisotropy 
in the meridional plane. The top-left panel shows the ra- 
dial over the angular velocity dispersion orIoq. This ra- 
tio does however not include the non-zero cross-term orq ■ 
The latter causes the velocity ellipsoid to be rotated with 
respect to the R and 9 coordinates. Taking this into ac- 
count the semi-axis lengths of the velocity ellipsoid in 
the meridional plane are given by a\ = (a 2 R + cr 2 )/2 ± 
yj (a R — <7g) 2 /4 + a Rg . In the top-right panel, we show the 
ratio of this minor <r_ and major cr + semi-axis length of 
the velocity ellipsoid (which is by definition in the range 
from zero to unity). This demonstrates that the velocity 
distribution of w Cen is nearly isotropic near the equato- 
rial plane, but becomes increasingly tangential anisotropic 
towards the symmetry axis. 

In the bottom panels we also include the azimuthal 
velocity dispersion o$. The bottom- left panel shows the 
radial over the tangential velocity dispersion, where the 
latter is defined as of = (of + o|)/2. Again this ratio 
does not take into account the cross-term grq. The ac- 
tual degree of anisotropy is given by the three semi-axis 
lengths o+, o_ and o^ of the velocity ellipsoid. In the 
bottom-right panel, we show, as a function of the position 
in the meridional plane, the minimum over the maximum 
of these three semi-axis lengths. Except for the region near 
the equatorial plane and within 10 arcmin, the best-fit 
model for lu Cen is clearly not isotropic. Even within this 
region, between about 3 and 5 arcmin, it is (slightly) ra- 
dially anisotropic. Outside this region lu Cen becomes in- 
creasingly tangentially anisotropic. 

Clearly, isotropic models are not suitable to model 
lu Cen. Also dynamical models with a two-integral dis- 
tribution function of the form F(E, L z ), with L z = R(v^) 
the angular momentum component along the symmetry z- 
axis, are not able to describe the complex dynamical struc- 
ture of lu Cen. For these models the solution of the Jeans 
equations can be used to construct dynamical models in 
a straightforward way (e.g. Satoh 1980; Binney, Davies & 
Illingworth 1990) and they allow for azimuthal anisotropy. 
However, for these models <jr — <jq and a R g = 0, i.e. 
isotropy in the full meridional plane, which is not the case 
for lu Cen (upper panels of Figure IT8"|) . Our axisymmetric 
dynamical models do not have these restrictions as they 
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Fig. 18. Degree of anisotropy as function of the equato- 
rial plane radius R and height z (excluding the axes to 
avoid numerical problems). The upper panels show the 
degree of anisotropy in the meridional plane: left the ra- 
dial over the angular velocity dispersion and right the mi- 
nor o_ over the major cr + semi-axis length of the veloc- 
ity ellipsoid, taking into account the cross-term a R g. The 
bottom panels include the azimuthal velocity dispersion: 
left the radial over the tangential velocity dispersion, with 
a 2 = (a 2 . + a 2 )/2, and right the minimum over the maxi- 
mum of the three semi-axis lengths o + , o_ and o^ of the 
velocity ellipsoid. See text for further details. 



are based on a general three-integral distribution func- 
tion F(E, L z , I3), which we investigate next for our best- 
fit model. 

9.3. Distribution function 

Each orbit in our models is characterised by the three in- 
tegrals of motion E, L z and I3. As function of these three 
integrals, we show in Figure 1191 for our best-fit model of 
lu Cen the distribution of the (mass) weights that were as- 
signed to the different orbits in the NNLS-fit. The energy 
E is sampled through the radius R c (in arcmin) of the cir- 
cular orbit (different panels) , of which we show the range 
that is constrained by the observations and that contains 
more than 90% of the total cluster mass. The angular mo- 
mentum L z (vertical) is in units of £ ma x, the angular mo- 
mentum of the circular orbit. The third integral Is (hor- 
izontal) is parameterised by the linearly sampled starting 
angle of the orbit, from the equatorial plane towards the 
symmetry axis, and of which the number is given. 

In each panel, the orbital weights are scaled with re- 
spect to the maximum orbital weight in that panel, indi- 
cated by the white colour, whereas black corresponds to 
zero orbital weight. The fraction of the sum of the mass 
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, 19. The orbital weight distribution for our best-fit model of a; Cen. From left to ri 
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Fig. 19. The orbital weight distribution for our best-fit model of a; Cen. From left to right, the panels show the orbital 
weight distribution at increasing distance from the centre, which corresponds to increasing energy. The radius R c (in 
arcmin) of the circular orbit at the corresponding energy is given above each panel. The radial range that is shown is 
constrained by the observations and contains more than 90% of the total cluster mass. The vertical axis represents the 
angular momentum L z in units of L maxi the angular momentum of the circular orbit. The horizontal axis represents 
the third integral I3, parameterised by the number of the (linearly sampled) starting angle of the orbit. Black shading 
corresponds to zero orbital weights, and white corresponds to the maximum orbital weight in each panel. At the 
bottom of each panel the fraction of the included mass with respect to the total mass is indicated (in %). 



weights in each panel with respect to total mass in all 
panels is given at the bottom of each panel (in %). To 
avoid an unrealistic orbital weight distribution that fluc- 
tuates rapidly for adjacent orbits, we regularise our mod- 
els fS l5.4|) . For values of the smoothening parameter below 
A = 4 and even without regularisation, we find the same 
best-fit parameters and although the distribution function 
becomes spiky, the main features of Figure HT^l remain. 

Most of the mass in the orbital weight distribution is 
in the component that is prominent in all panels. With 
increasing radius, the average angular momentum L z of 
this component increases from nearly zero to a significant 
(positive) value in the outer parts. This reflects the out- 
wards increasing tangential anisotropy already seen in the 
bottom-left panel of Figure ED An almost non-rotating 
part is still present beyond 5 arcmin, attached to the ro- 
tating component, which becomes the dominant compo- 
nent (in mass). There is also a separate component at 
L z /L max ~ 1 that is clearly visible between about 1 and 
3 arcmin. Within this radial range, this maximum rotat- 
ing component contributes almost 20% of the mass, and 
it includes about 4% of the total mass, i.e., its mass is of 
the order of 10 5 M Q . 

In the right-most panels of Figure HIH there is a (weak) 
signature of a component with i z /L max ~ —1, which we 
expect to be a spurious feature due to insufficient ob- 
servational constraints. Whereas (nearly) circular orbits 
(|_L Z |/L max ~ 1) are confined in radius to R c , orbits with 
lower \L Z \ can go further inwards, so that they have most 
of their contribution (their cusps) at a smaller radius than 
R c (e.g. Cappellari et al. 2004). Hence, the apparent fea- 
ture at L z /L max ~ — 1 in the most-right panel is only con- 
strained by data around and beyond the radius R c — 13.6 



arcmin, where the coverage of the data is sparse with only 
a few polar bins (see Figure 113(1 . The main component 
in this panel at L z /L max ~ 0.5 is (mostly) constrained 
by data at smaller radii, where there is good data cover- 
age. The separate maximum rotating component between 
1 and 3 arcmin is constrained by only a few proper motion 
apertures, but is strongly constrained by the line-of-sight 
velocity data. 

Due to the difference in spatial coverage between the 
proper motion and line-of-sight velocity data, the two 
data-sets (better) constrain different parts of the orbital 
weight distribution. By fitting besides the light distribu- 
tion of w Cen the mean velocity and velocity dispersion of 
only the proper motion components, we find a less promi- 
nent separate component between 1 and 3 arcmin, but it 
is still present. In the case of only fitting the mean line-of- 
sight velocity and velocity dispersion, this separate com- 
ponent is clearly visible and even extends into the outer ro- 
tating main component. The transition between the main 
non-rotating and rotating component is in the case of only 
line-of-sight data more abrupt than in Figure^J However, 
the proper motion data, which has a better coverage in the 
outer parts, shows a similar smooth transition. We con- 
clude that, although the spatial coverage is different, both 
data-sets give rise to the same main features in the orbital 
weight distribution. 

9.4. Dynamical substructures 

Within 5 arcmin the main component has on average a 
high value of I3. In combination with the low value of L z , 
we interpret this as a non-rotating spheroidal structure. 
Beyond 5 arcmin, L z increases and I3 decreases, and the 
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main component flattens and rotates faster. The smaller 
component attached to it may well be the signature of the 
fading non-rotating spheroidal component. 

For the separate component between 1 and 3 arcmin, 
L z approaches its maximum value. As a result, the zero- 
velocity curve shrinks towards the circular orbit in the 
equatorial plane, and the corresponding orbits are all flat, 
irrespective of the (high) value of 23 (see also Figure 3 of 
Cretton et al. 1999). Hence, this fast-rotating component 
is likely to be an inner disk, which fades away into the 
more massive main rotating component at larger radius. 

We compute the spatial distribution and average kine- 
matics of these possible substructures in the phase-space 
of uj Cen. To this end we select the orbits from our best-fit 
model that contribute non-zero weight to three different 
parts of the distribution function in Figure We select 
the inner main component in the 7 left-most panels, ex- 
cluding the separate disk component in the 5 left-most 
panels, and the outer main component in the 3 right- 
most panels (excluding the weak feature in the bottom). 
For each orbit with non-zero weight, we then randomly 
draw points along its numerically integrated orbit, with 
the number of drawings proportional to its relative weight. 
In this way, we make an (N-body) realisation of our best-fit 
model consisting of a couple of tens of thousands of par- 
ticles, representing the stars in u> Cen. For each of these 
stars, we determine the position on the plane of the sky 
and the three velocity components; the two proper motion 
components in the plane of the sky and the line-of-sight 
velocity. For the stars that belong to a certain part or 
substructure of phase-space, we then calculate the spatial 
distribution and mean kinematics. 

Figure I2U1 shows the results for all stars, those in the 
inner and outer main component and those in the sepa- 
rate disk component, respectively, per column from left 
to right. The first row shows the spatial distribution. The 
flattening of the spatial distribution of all stars and of the 
outer main component are both about 0.88, similar to the 
average observed flattening for uj Cen. The inner main 
component, going out to a radius of about 6 arcmin, is 
rounder with a flattening of about 0.94. The spatial dis- 
tribution of the disk component only extends to a radius of 
about 3 arcmin, has an average flattening as lows as 0.60 
and is less dense in the centre as this maximum rotat- 
ing disk consists of stars on (nearly) circular orbits which 
avoid the centre. The second to fourth row show the mean 
velocity fields in respectively the direction of the major x ' - 
axis and the minor y'-axis on the plane of the sky and the 
line-of-sight z'-axis. In each panel the axes are scaled with 
respect to the spatial extent of each component. Whereas 
the inner main component indeed hardly show any rota- 
tion, the outer main component clearly rotates and the 
separate disk component rotates even faster. In the last 
row, the velocity dispersion profiles are presented, radial 
(green) and tangential (red) on the plane of the sky and 
along the line-of-sight (blue) . Even though the outer main 
component is flatter and rotates faster than the inner main 
component, it is not kinematically colder due to the mix- 



ture of orbits with different h z values. On the other hand, 
the maximum rotating disk is the kinematically coldest 
component. Whereas the inner main component is nearly 
isotropic, the outer main component is anisotropic and the 
disk component is even stronger anisotropic. 

The presence of dynamical substructures implies that 
the formation history of uj Cen is more complicated than 
expected for a typical globular cluster. However, the inter- 
pretation of these different components in the distribution 
function is very difficult. In what follows we investigate the 
possible effects due to the tidal interaction between uj Cen 
and the Milky Way I9.5JI . and the possible link to the 
observed multiple stellar populations in uj Cen fg 19.6(1 . 

9.5. Tidal interaction 

Based on its current position and motion in the Milky 
Way (MW), Dinescu et al. (1999) simulated the orbit of 
uj Cen around the Galactic Centre (GC). They found that 
the average orbit is inclined by only 17° with respect to 
the Galactic plane, has a period of P ~ 122 Myr and an 
angular momentum of about 406 kpc km s" 1 . Assuming 
that the average orbit of uj Cen is circular, we thus find a 
radius Roc ~ 2.8 kpc and a velocity of about 143 kms" 1 , 
of which the component perpendicular to Galactic plane 
v± ~ 42 kms -1 . Since the scale height of the MW disk is 
typically 250 pc, it takes about t cnc ~ 12 Myr for uj Cen 
to cross the MW disk. This means that for nearly 10% 
of its time uj Cen is immersed in the disk and feels the 
additional gravitational field. 

Based on its current position and motion in the Milky 
Way (MW), Dinescu et al. (1999) simulated the orbit of 
uj Cen around the Galactic Centre (GC). They found that 
the average orbit is eccentric (e ~ 0.68), is inclined by 
only 17° with respect to the Galactic plane and has a 
period of P ~ 122 Myr. In combination with the orbital 
angular momentum of about 406 kpc km s" 1 , this implies a 
mean circular radius of Roc ~ 2.8 kpc and corresponding 
circular velocity of 143 km s _1 . The velocity perpendicular 
to Galactic plane is on average thus v± ~ 42 kins" 1 . Since 
the scale height of the MW disk is typically 250 pc, it takes 
about t cnc ~ 12 Myr for uj Cen to cross the MW disk. This 
means that for nearly 10% of its time uj Cen is immersed 
in the disk and feels the additional gravitational field. 

To investigate what effect the MW tidal field has on 
the stars in uj Cen, we use the impulse approximation as 
described by Binney & Tremaine (1987, p. 446), with the 
typical properties of the MW from their Tables 1-1 and 1- 
2. We assume a Cartesian coordinate system with its ori- 
gin at the centre of uj Cen and the z-axis perpendicular 
to the MW disk. If uj Cen goes through the MW disk, the 
effect on the velocity component perpendicular to the disk 
is the largest. Hence, the velocity of a cluster star changes 
on average by |Au| ~ z\g z (R)\/v±, where g z is the z- 
component of the gravitational field of the MW disk. The 
cumulative effect of successive passages through the MW 
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Fig. 21. Timescales as function of the projected radius 
R' . The green curve represents the timescale on which 
shocks, caused by successive passages of to Cen through 
the MW disk, change the mean-squared velocity of a clus- 
ter star by the order of the (local) velocity dispersion of 
the cluster. The red and blue curves show respectively the 
dynamical time tdyn and relaxation time £ re iax- The verti- 
cal dashed lines indicate with increasing distance the core 
radius r*, the half-light radius Th and the tidal radius rt 
of ui Cen. 

disk becomes of the order of the (local) velocity dispersion 
g of the cluster on a timescale of i s hock ~ Po 2 v\j (8z 2 g 2 ). 

An infinite disk with surface density £ generates a 
gravitational field g z — 2-kGYj. In the solar neighbourhood 
the MW disk has a surface density of E ~ 75 M Q pc~ 2 . 
Assuming that the MW disk falls off as exp(—R/Rd) in the 
radial coordinate, with Rd = 3.5 kpc , we find that at the 
mean circular radius R = Roc of to Cen's orbit around 
the GC, g z ~ 2.9 x 10~ 13 kms~ 2 . For a spherical shell 
of stars of radius r, we have that on average z 2 = r 2 /3. 
We thus find that the timescale on which disk shocking 
becomes important is 




Figure [21 shows £ s hock (green curve) as function of the 
projected radius R' (in arcmin). We used the line-of-sight 
velocity dispersion as given in Figure |H1 smoothed and 
extrapolated to larger radii with the help of measurements 
by Scarpa, Marconi & Gilmozzi (2003) between about 20 
and 30 arcmin 8 . In the same figure we have also plotted 

8 Taking into account the measurement error of about 1 
kms -1 and the perspective rotation that can be as large as 
1.5 kms -1 at those radii (eq. |SJ. 



the dynamical time idyn (red curve; Binney & Tremaine 
1987, eq. 2-30) and the relaxation time i re iax (blue curve; 
Spitzer & Hart 1971; Binney & Tremaine 1987, eq. 8-71). 
The three vertical dashed lines indicate respectively the 
core radius r c = 2.6 arcmin, the half-light radius 77, = 4.8 
arcmin and the the tidal radius r t = 45 arcmin (e.g. Trager 
et al. 1995). 

Clearly, the impulse approximation is not valid near 
the centre of ui Cen, where the period of the stellar or- 
bits T = 4t c iyn is much smaller than the duration of the 
passage through the disk t cnc ~ 12 Myr. Disk shocking is 
thus unimportant at the centre of u> Cen: the orbits evolve 
adiabatically and emerge unharmed from the encounter. 
Around a radius of 16 arcmin, where T is about twice t enc , 
disk shocks begin to play an important role since the disk 
shocking time becomes of the order of the dynamical time 
£shock ~ tdyn ~ 6 Myr. At the tidal radius of 45 arcmin, 
the MW disk gravitational field becomes dominant. 

The effect that the MW tidal field has on the internal 
dynamics of u> Cen also strongly depends on the relative 
orientation and spinning direction of the angular momen- 
tum vector of the stars in ui Cen (internal) and the angu- 
lar momentum vector of its orbit around the GC (exter- 
nal). We found that the rotation axis is about 50° inclined 
with respect the line-of-sight (the z'-axis) in the direction 
South 9 . On the plane of sky, the rotation axis projects onto 
the minor j/'-axis, which makes an angle of about 10° away 
from North in the direction East. The equatorial coordi- 
nates of lu Cen are a = 13 ft 26 m 46 s and S = -47°28'43" 
(J2000), which correspond to a Galactic longitude and lat- 
itude of I = 309° and b — 15°. Hence, the rotation axis is 
nearly parallel (angle < 3°) to the equatorial plane, and 
makes an angle of about 65° with respect to the Galactic 
plane. Seen from the North Galactic pole, lu Cen is moving 
in anti-clockwise direction around the GC. The rotation 
inside lu Cen is dominated by orbits with positive L z val- 
ues in Figure ^3 which correspond to clockwise rotation. 

We thus find that the internal and external angular 
momentum vector are for more than 90% parallel with re- 
spect to each other with opposite spinning direction. From 
mergers of spinning galaxies it is well known that if the 
spins are anti-parallel as in this case, the orbital disrup- 
tion is much less than in the case of parallel spins (e.g. 
Toomre & Toomre 1972). Hence, in the past lu Cen might 
have contained a significant number of stars on orbits with 
negative L z (parallel spin) , which then were removed from 
the cluster during its successive passages through the MW 
disk. On the other hand, stars on orbits with positive L z 
(anti-parallel spin) had a bigger chance to survive. 

Furthermore, the stars on more radial orbits (those 
with smaller values of L z ) cover a broader range in radius, 



This means that in the common definition of the inclina- 
tion, as in eq. the best-fit inclination is -50°. This also 
explains the sign difference of (tv) in eq. @ and along the ver- 
tical axis of the plot in the middle panel of Figure [7] However, 
we decided to adopt the usual convention to take the value for 
the inclination in the range from 0° (face-on) to 90° (edge-on). 
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with the influence of the MW tidal field becoming stronger 
at increasing radius. In the course of time, these radial 
orbits thus have a bigger chance of being disrupted than 
the more tangential orbits with similar mean radius. 

Both effects (together) might explain the prominent 
rotating main component in the distribution function in 
Figure 1191 beyond a radius of 10 arcmin, while the non- 
rotating main component that dominates inwards, fades 
away. The removal of the more radial orbits also naturally 
explains the outwards increasing tangential anisotropy in 
our best-fit model of w Cen (§ 19. 2|) . 

The above analysis shows that the frequent passages 
of lo Cen through the MW disk most likely have played a 
crucial role in the evolution of this cluster. At least part of 
the phase-space structure of lo Cen may well be caused by 
the tidal field of the MW. Detailed (N-body) simulations 
are needed to further quantify this. 

9.6. Multiple stellar populations 

Among the Galactic globular clusters, lo Cen especially 
stands out because of its chemical inhomogeneity, first re- 
vealed in photometric investigations by Dickens & Woolley 
(1967) and spectroscopically confirmed by Freeman & 
Rodgers (1975). Besides the main population of metal- 
poor stars (~65% of all stars with [Ca/H]~ —1.4) and an 
intermediate population (~30%, [Ca/H]~ —1.0), recently 
also a separate metal-rich population (~5%, [Ca/H]~ 
—0.5) has been identified (Lee et al. 1999; Pancino et al. 
2000), and even the main sequence of lo Cen is bifurcated 
(Bedin et al. 2004). 

Theses different stellar populations also appear to have 
a different spatial distribution. Whereas the metal-poor 
stars seems to follow the observed flattening of lo Cen 
in the East- West direction, the more metal-rich stars are 
elongated in the North-South direction and also more cen- 
trally concentrated (e.g. Pancino et al. 2003). There are 
also indications of differences in the kinematics of the stel- 
lar populations. Norris et al. (1997) find that the metal- 
poor populations have on average a higher line-of-sight 
velocity dispersion and exhibit a well-defined linc-of-sight 
rotation, while the metal-rich populations show no signif- 
icant rotation. Ferraro, Bellazzini & Pancino (2002) claim 
that the separate metal-rich population has a coherent 
bulk proper motion significantly different from the other 
cluster stars. 

We use the empirical relation in eq. (15) of Paper I to 
estimate the [Ca/H] abundances of stars in our analysis 
with l/-band magnitude and B — V colour measurements 
consistent with the top of the red giant branch (V < 13.5 
and B — V > 0.7). The resulting [Ca/H] histograms for the 
proper motion and line-of-sight velocity stars both show 
a distribution with a broad peak around [Ca/H]^ —1.2 
and a long tail extending beyond [Ca/H]^ —0.5. In both 
cases the peak shows a small dip, so that we might divide 
the stars into a metal-poor population with [Ca/H]< —1.2 



and a metal-rich population with [Ca/H]> —1.2, similar 
to Norris et al. (1997). 

Comparing the mean line-of-sight kinematics of the 
metal-poor and metal-rich stars, we confirm the result of 
Norris et al. (1997) that the more centrally concentrated 
metal-rich stars are on average kinematically cooler and 
nearly non-rotating. The line-of-sight velocity dispersion 
profile is steeper for the metal-richer stars than for the 
metal-poor stars, such that that in the centre the metal- 
richer stars are even (slightly) kinematically warmer. The 
proper motions seems to imply a similar difference in the 
slope of the velocity dispersion profiles. However, with the 
proper motion errors on average four times larger than 
those of the line-of-sight velocities (see also Figure |HJ, 
there are no significant differences between the kinematics 
of the metal-poor and metal-rich stellar populations. 

The above correlations between the kinematics and 
chemical properties of stars in lo Cen, are expected to 
show up in the distribution function (see also Freeman 
2002). The centrally concentrated non-rotating metal- rich 
stars would lie near the bottom of the potential well at 
the lower values of E found in the cluster, symmetrically 
distributed over positive and negative values of L Z: and 
towards higher values of I3 . The rotating metal-poor stars 
would span the entire range of E, with an asymmetric 
distribution in L z and towards lower I3. 

These expectations are consistent with the orbital 
weight distribution of our best-fit dynamical model of 
lo Cen fFigure Hill and l2"0|l . Whereas the metal-richer stars 
might well be associated with the inner non-rotating part 
of the main component, we might see the kinematical 
signatures of the metal-poorer stars becoming dominant 
when the main component flattens and rotates faster in 
the outer parts. Still, we have to be careful as these are (in- 
direct) indications of a link between substructures in the 
distribution function and the different stellar populations. 

To investigate directly the distribution function of the 
different stellar populations, once can try to construct sep- 
arate dynamical (Schwarzschild) models. However, since 
the separation into different stellar populations is not ev- 
ident, separate mass models are needed and the separate 
kinematic constraints are based on much fewer stars, this 
is very difficult with the current data-set. A more feasi- 
ble approach would be to model together, in a consistent 
way, the observed kinematics and physical properties of 
the stars. For example, by labelling the orbits in the model 
with different colours, the observed colour (averaged per 
aperture) can be used to constrain the model in addition 
to the photometry and kinematics. On the other hand, 
now that we have constrained the global parameters (dis- 
tance, inclination and mass-to-light ratio) considerably, it 
has become feasible to use non-linear maximum likelihood 
techniques to directly incorporate discrete stellar measure- 
ments. In this way, for the model that best fits (simulta- 
neously) the measured kinematics and age and metallicity 
indicators of individual stars, the different stellar popula- 
tions can be cleanly separated in phase-space. This exten- 
sion, which we leave for a future paper, will provide an 
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important contribution to solving the stellar population 
puzzle in lo Cen, and clarify its formation history. 

10. Conclusions 

We used an extension of Schwarzschild's (1979) orbit su- 
perposition method to construct realistic axisymmetric 
dynamical models for lo Cen with an arbitrary anisotropic 
velocity distribution. By fitting these models simultane- 
ously to proper motion and line-of-sight velocity measure- 
ments, we measured the radial mass-to-light profile, the 
inclination and the distance to lo Cen, which is needed to 
convert the proper motions to physical units. This dynam- 
ical distance estimate can provide a useful calibration for 
the photometric distance ladder. 

We used the ground-based proper motions from 
Paper I and the line-of-sight velocities from four indepen- 
dent data-sets. We brought the kinematic measurements 
onto a common coordinate system and carefully selected 
on cluster membership and on measurement error. This 
provided a homogeneous data-set of 2295 stars with proper 
motions accurate to 0.20 masyr -1 and 2163 stars with 
line-of-sight velocities accurate to 2 kms -1 , covering a ra- 
dial range out to about half the tidal radius of the cluster. 
We corrected the kinematic measurements for perspective 
rotation and removed a residual solid-body rotation com- 
ponent in the proper motions. We showed that the latter 
can be measured without any modelling other than assum- 
ing axisymmetry and at the same time we obtained a tight 
constraint on Dtani of 5.6 (+1.9/-1.0) kpc, providing a 
unique way to estimate the inclination i of a nearly spher- 
ical object once the distance D is known. The corrected 
mean velocity fields are consistent with regular rotation, 
and the mean velocity dispersions display significant de- 
viations from isotropy. 

We binned the individual measurements on the plane 
of the sky to search efficiently through the parameter space 
of the models. Tests on an analytic model demonstrated 
that our approach is capable of measuring the cluster dis- 
tance to an accuracy of about 6 per cent. Application 
to lo Cen revealed no dynamical evidence for a signifi- 
cant radial dependence of the (V-band) stellar mass-to- 
light ratio M/Ly, in harmony with the relatively long 
relaxation time of the cluster. We found that our best- 
fit dynamical model has M/Ly = 2.5 ± 0.1 M©/L© and 
i = 50° ± 4°, which corresponds to an average intrinsic 
axial ratio of 0.78 ± 0.03. The best-fit dynamical distance 
D = 4.8 ± 0.3 kpc (distance modulus 13.75 ± 0.13 mag) 
is significantly larger than obtained by means of simple 
spherical or constant-anisotropy axisymmetric dynamical 
models, and is consistent with the canonical value 5.0±0.2 
kpc obtained by photometric methods. The total mass of 
the cluster is (2.5 ± 0.3) x 10 6 M Q . 

Schwarzschild's approach also provides an insight into 
the intrinsic orbital structure of the cluster. Our best-fit 
model implies that lo Cen is close to isotropic inside a 
radius of about 10 arcmin and becomes increasingly tan- 
gentially anisotropic in the outer region, which displays 



significant mean rotation. We found that this may well 
be caused by the effects of the tidal field of the Milky 
Way. Furthermore, the best-fit model contains a separate 
disk-like component between 1 and 3 arcmin, contributing 
about 4% to the total mass. This phase-space structure, 
which might be linked to the multiple stellar populations 
observed in lo Cen, is expected to provide important con- 
straints on its formation history. 

We might improve our best-fit dynamical model of 
lo Cen and better constrain the distance and the other 
parameters, by extending the data-set with for example 
proper motions derived from HST images. Whereas with 
the ground-based proper motions we were unable to probe 
the centre of lo Cen due to crowding, the high spatial reso- 
lution and high sensitivity of HST, results in many proper 
motion measurements in the very centre, which makes it 
possible to investigate a possible central mass concentra- 
tion in lo Cen. 

We may also increase the kinematic constraints on 
our dynamical models by including mean correlated and 
higher-order velocity moments. With the parameter range 
considerably constrained, it now becomes also feasible to 
use non-linear maximum likelihood techniques to directly 
incorporate the discrete kinematic measurements. These 
techniques not only allow correlated and higher-order ve- 
locity moments to be included in a straightforward way, 
but also provide a natural way to incorporate measure- 
ments of age and metallicity indicators of individual stars 
in addition to their photometry and kinematics. By fit- 
ting an orbit-based model simultaneously to all these ob- 
servations, different stellar populations can be separated 
in phase-space, after which their structure and dynamics 
can be studied separately. 

We have shown that with the method described in this 
paper, we were able to measure the global parameters of 
lo Cen, including its distance, and investigate its intrin- 
sic orbital structure. This method can also be applied to 
study other globular clusters and stellar clusters in the 
Milky Way, provided that accurate velocity measurements 
are available. With the amount of (photometric and kine- 
matic) data quickly increasing, we expect this method to 
become an important tool to model these stellar systems 
and gain insight in their formation and evolution. 

Acknowledgements 

This project would not have been possible without the 
extensive plate collection of lo Cen obtained in the nine- 
teen thirties by Willcm Christiaan Martin, which, together 
with the plates taken in the eighties by Ken Freeman and 
Pat Seitzer and the massive subsequent effort by Rudolph 
Le Poole and Floor van Leeuwen, produced the proper 
motions reported in Paper I. The line-of-sight velocities 
used here include those obtained at CTIO in the early 
nineties by Renate Reijns and Pat Seitzer, as reported in 
Paper II. We are also grateful to Karl Gcbhardt for allow- 
ing us to use his unpublished Fabry-Perot measurements, 
and to Jay Anderson & Ivan King for providing their pre- 



32 



G. van de Ven et al.: The dynamical distance and intrinsic structure of u> Cen 



liminary HST proper motions. It is a pleasure to thank 
Thijs Kouwenhoven for assistance in the initial phase of 
the work reported here, to thank Jesus Falcon-Barroso, 
Ken Freeman, Floor van Leeuwen, and Scott Tremaine 
for useful discussions and suggestions during the course 
of this work, and to thank Michele Cappellari, Richard 
McDermid and Anne-Marie Weijmans for a critical read- 
ing of the Omanuscript. We also thank the referee for con- 
structive comments and suggestions that improved the 
presentation of the paper. This research was supported 
in part by NWO through grant 614.000.301, and by the 
Netherlands Research School for Astronomy NOVA. 

References 

Bedin, L. R., Piotto, G., Anderson, J., et al. 2004, ApJ, 605, 
L125 

Benedict, G. F., McArthur, B. E., Fredrick, L. W., et al. 2002, 
AJ, 123, 473 

Binney, J. & Tremaine, S. 1987, Galactic Dynamics (Princeton, 

NJ, Princeton University Press) 
Binney, J. J., Davies, R. L., & Illingworth, G. D. 1990, ApJ, 

361, 78 

Cappellari, M. 2002, MNRAS, 333, 400 

Cappellari, M., van den Bosch, R. C. E., Verolme, E. K., et al. 
2004, in Carnegie Observatories Astrophysics Series, Vol. 
1, Coevolution of Black Holes and Galaxies, ed. L. C. Ho 



Gerhard, O. E. 1993, MNRAS, 265, 213 

Geyer, E. H., Nelles, B., & Hopp, U. 1983, A&A, 125, 359 

Harris, W. E. 1996, AJ, 112, 1487 

Hog, E., Fabricius, C, Makarov, V. V., et al. 2000, VizieR 

Online Data Catalog, 1259, 
Icke, V. & Alcaino, G. 1988, A&A, 204, 115 
Kalnajs, A. J. 1999, in ASP Conf. Ser. 165: The Third Stromlo 

Symposium: The Galactic Halo, 325 
Kenney, J. F. & Keeping, E. S. 1951, The Distribution of the 

Standard Deviation (§ 7.8 in Mathematics of Statistics, Pt. 

2, 2nd ed. Princeton, NJ: Van Nostrand) 
King, I. R. & Anderson, J. 2002, in ASP Conf. Ser. 265: Omega 

Centauri, A Unique Window into Astrophysics, eds. F. van 

Leeuwen, J. D. Hughes, G. Piotto, 21 
King, I. R., Hedemann, E. J., Hodge, S. M., & White, R. E. 

1968, AJ, 73, 456 
Kleyna, J., Wilkinson, M. I., Evans, N. W., Gilmore, G., & 

Frayn, C. 2002, MNRAS, 330, 792 
Krajnovic, D., Cappellari, M., Emsellem, E., McDermid, R. M., 

& de Zeeuw, P. T. 2005, MNRAS, in press 
Lawson, C. L. & Hanson, R. J. 1974, Solving least squares 

problems (Prentice-Hall Series in Automatic Computation, 

Englewood Cliffs: Prentice-Hall, 1974) 
Lee, Y.-W., Joo, J.-M., Sohn, Y.-J., et al. 1999, Nature, 402, 

55 

Lee, Y.-W., Rey, S.-C, Ree, C. H., et al. 2002, in ASP 
Conf. Ser. 265: Omega Centauri, A Unique Window into 
Astrophysics, eds. F. van Leeuwen, J. D. Hughes, G. Piotto, 
3Q5 



(http^/www.ociw.edu/ociw/symposia/series/symposiuml/pr^edipg^J^ml^ Agp Conf ger 265 . 0mega Centauri5 A 
Cappellari, M., Verolme, E. K, van der Marel, R. P., et al. 



Rix, H. 



2002, ApJ, 578, 787 
Carraro, G. & Lia, C. 2000, A&A, 357, 977 
Cretton, N., de Zeeuw, P. T., van der Marel, R. 1 

1999, ApJS, 124, 383 
Cudworth, K. M., 1979, AJ, 84, 1312 
Da Costa, G. S. 1979, AJ, 84, 505 

Davies, R. L., Efstathiou, G., Fall, S. M., Illingworth, G 

Schechter, P. L. 1983, ApJ, 266, 41 
Dejonghe, H. 1989, ApJ, 343, 113 

Dickens, R. J. & Woolley, R. v. d. R. 1967, Royal Greenwich 

Observatory Bulletin, 128, 255 
Dinescu, D. I., Girard, T. M., & van Altena, W. F. 1999, AJ, 

117, 1792 

Dull, J. D., Cohn, H. N., Lugger, P. M., et al. 1997, ApJ, 481, 
267 

Emsellem, E., Monnet, G., & Bacon, R. 1994, A&A, 285, 723 
Emsellem, E., Monnet, G., Bacon, R., & Nieto, J.-L. 1994, 

A&A, 285, 739 
Evans, N. W. 1994, MNRAS, 267, 333 
Evans, N. W. & de Zeeuw, P. T. 1994, MNRAS, 271 
Feast, M. W., Thackeray, A. D., & Wesselink, A. 

MNRAS, 122, 433 
Ferraro, F. R., Bellazzini, M., & Pancino, E. 2002, ApJ, 573, 

L95 

Freeman, K. C. 1993, in Astronomical Society of the Pacific 

Conference Series, 608 
Freeman, K. C. 2002, in ASP Conf. Ser. 265: Omega Centauri, 

A Unique Window into Astrophysics, eds. F. van Leeuwen, 

J. D. Hughes, G. Piotto, 423 
Freeman, K. C. & Rodgers, A. W. 1975, ApJ, 201, L71 
Gascoigne, S. C. B. & Burr, E. J. 1956, MNRAS, 116, 570 
Gebhardt, K, Richstone, D., Tremaine, S., et al. 2003, ApJ, 

583, 92 



202 

J. 1961, 



Unique Window into Astrophysics, eds. F. van Leeuwen, 
J. D. Hughes, G. Piotto, 95 
Lucy, L. B. 1974, AJ, 79, 745 

Mandushev, G., Staneva, A., & Spasova, N. 1991, A&A, 252, 
94 

Martin, W. C. 1938, Annalen van de Sterrewacht te Leiden, 
Vol. XVII 

Mayor, M., Duquennoy, A., Udry, S. et al. 1996, in ASP Conf. 
Ser. 90: The Origins, Evolution, and Destinies of Binary 
Stars in Clusters, eds. E. F. Milone, J.-C. Mermilliod, 190 

Mayor, M., Meylan, G. , Udry, S., et al. 1997, AJ, 114, 1087 
[M97] 

McNamara, B. J., Harrison, T. E., Baumgardt, H., 2004, ApJ, 
602, 264 

Merrifield, M. R. & Kent, S. M. 1990, AJ, 99, 1548 
Merritt, D. 1993, ApJ, 413, 79 
— . 1997, AJ, 114, 228 

Merritt, D., Meylan, G., & Mayor, M. 1997, AJ, 114, 1074 
Merritt, D. & Saha, P. 1993, ApJ, 409, 75 
Meylan, G. 1987, A&A, 184, 144 

Meylan, G., Mayor, M., Duquennoy, A., & Dubath, P. 1995, 
A&A, 303, 761 

Monnet, G., Bacon, R., & Emsellem, E. 1992, A&A, 253, 366 
Norris, J. E., Freeman, K. C, Mayor, M., & Seitzer, P. 1997, 
ApJ, 487, L187 

Pancino, E., Ferraro, F. R., Bellazzini, M., Piotto, G., & 

Zoccali, M. 2000, ApJ, 534, L83 
Pancino, E., Seleznev, A., Ferraro, F. R., Bellazzini, M., & 

Piotto, G. 2003, MNRAS, 345, 683 
Perryman, M. A. C. & ESA. 1997, The HIPPARCOS and 

TYCHO catalogues (ESA SP 1200) 
Peterson, R. C, Cudworth, K. M., 1994, ApJ, 420, 612 
Peterson, R. C, Rees, R. F., Cudworth, K. M., 1995, ApJ, 443, 

124 



G. van de Ven et al.: The dynamical distance and intrinsic structure of u) Cen 



33 



Platais, I., Wyse, R. F. G., Hebb, L., Lee, Y., & Rey, S. 2003, 
ApJ, 591, L127 

Press, W. H., Teukolsky, S. A., Vettering, W. T., & Flannery, 
B. P. 1992, Numerical Recipes (Cambridge Univ. Press, 
Cambridge) 

Rees, R. F., 1997, in ASP Conf. Ser. 127: Proper Motions and 
Galactic Astronomy, ed. R.M. Humphreys, 109 

Reijns, R. A., Seitzer, P., Arnold, R., et al. 2005, A&A, sub- 
mitted [Paper II] 

Richstone, D. O. & Tremaine, S. 1988, ApJ, 327, 82 

Rix, H., de Zeeuw, P. T., Cretton, N., van der Marel, R. P., & 
Carollo, C. M. 1997, ApJ, 488, 702 

Romanowsky, A. J. & Kochanek, C. S. 2001, ApJ, 553, 722 

Satoh, C. 1980, PASJ, 32, 41 

Scarpa, R., Marconi, G., & Gilmozzi, R. 2003, A&A, 405, L15 
Schwarzschild, M. 1979, ApJ, 232, 236 
— . 1982, ApJ, 263, 599 
— . 1993, ApJ, 409, 563 

Seitzer, P. O. 1983, Ph.D. Thesis, Australian National Univ. 
Spitzer, L. J. & Hart, M. H. 1971, ApJ, 164, 399 
Suntzeff, N. B. & Kraft, R. P. 1996, AJ, 111, 1913 [SK96] 
Thompson, I. B., Kaluzny, J., Pych, W., et al. 2001, AJ, 121, 
3089 

i Toomre, J. 1972, ApJ, 178, 623 
, King, I. R., & Djorgovski, S. 1995, AJ, 109, 218 
Gebhardt, K., Bender, R., et al. 2002, ApJ, 574, 



Toomre, A. 
Trager, S. C 
Tremaine, S. 
740 

Tsuchiya, T., Korchagin, V. L, & Dinescu, D. I. 2004, MNRAS, 
350, 1141 

van der Marel, R. P. 2004, in Carnegie Observatories 
Astrophysics Series, Vol. 1, Coevolution of 
Black Holes and Galaxies, ed. L. C. Ho 
l |http://www.ociw.edu7o ciw/symposia/series/symposiuml/proceedings. n"tml^ 



A possible way to measure the mean velocity and ve- 
locity dispersion, is to fit a Gaussian distribution to the 
velocity histogram of the stars that fall within an aperture. 
Whereas the mean velocity V is well estimated, the best-fit 
mean velocity dispersion am is too large, as the Gaussian 
distribution is broadened due to the velocity errors. This 
additional 'instrumental' dispersion cri ns can be estimated 
by the mean of the velocity errors. The corrected mean 
velocity dispersion a then follows from a 2 = <j\ t — af ns . 
Since this is only an approximate correction, we use a max- 
imum likelihood estimate of the velocity moments that at 
the same time corrects for each individual velocity error. 

Suppose L{v ) is the (intrinsic) velocity distribution of 
the stars in an aperture, in one of the three principal di- 
rections. We can consider each stellar velocity measure- 
ment Vi in that aperture as drawn from this distribution, 
or alternatively, the product of L(v) with a delta func- 
tion around Uj, integrated over all velocities. Due to (in- 
strumental) uncertainties this delta-function is broadened, 
and we assume that it can be described by a Gaussian 
around Vi, with the corresponding velocity error Cj as the 
standard deviation. For a sufficient number of draws N, 
i.e. velocity measurements in the aperture, we can then 
recover the (unknown) velocity distribution C(v) by max- 
imising the likelihood 
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Appendix A: Maximum likelihood estimation 
velocity moments 

We use the average kinematics of stars that fall within 
apertures on the plane of the sky. This is comparable to 
the kinematics from the integrated spectra of galaxies in 
an aperture. A very important difference is, however, that 
we have to take into account the errors on the individual 
velocity measurements. 



or, equivalently, minimising A = — 21n(L), with respect 
to the mean velocity V, mean velocity dispersion a and 
possible higher-order velocity moments. 

It is possible to recover C(v) in a non-parametric way 
using (extensions of) Lucy's (1974) method, but exploit- 
ing the fact that Gaussians are good low-order approxi- 
mations, the velocity distribution is often parameterised 
by a Gauss-Hermite (GH) scries (van der Marel & Franx, 
1993; Gerhard, 1993). It has the advantage that it only re- 
quires the storage of the velocity moments (V, a, ft.4, 
. . . ) instead of the full velocity distribution. Furthermore, 
it allows a simple velocity scaling of the model, which is 
useful when investigating the effect of a change in the stel- 
lar mass-to-light ratio. 

Another advantage of parameterising C(v) comes from 
the observation that the integral in IjA.ip is the convolu- 
tion of the velocity distribution and the Gaussian of each 
velocity measurement. For a Gaussian velocity distribu- 
tion this convolution is straightforward, but also in the 
case that C(v) is described by a GH series, the convolution 
can be carried out analytically. This makes it feasible to 
apply the method to a large number of discrete measure- 
ments and to estimate the uncertainties on the extracted 
velocity moments by means of the Monte Carlo bootstrap 
method (§ 15.6 of Press et al. 1992). 

In the case of no measurement errors, the maximum 
likelihood estimator of the standard deviation a, given by 
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proper motions 



\ 



1 

n ' 

i=i 



with v = y^Vi, (A. 2) 

i=l 



is a biased estimator, underestimating the true a by a 
factor (see also e.g. Kenney & Keeping 1951, p. 171) 



b{n) 



2 r(f) 



4n 32n 2 



(A.3) 



where T is the gamma function. When we take into ac- 
count the measurement errors <7j, there is no such simple 
analytical bias correction as (|A.3|I . However, we can use 
the latter result to derive the following approximate cor- 
rected standard deviation estimator 



1 



b(n) 



<r 2 + [l-b 2 (n)]a 2 , 



(A.4) 



where a is the maximum likelihood estimated dispersion 
and a 2 = ^ $^i=i °f ^ ne average measurement error. 

Appendix B: Polar grid of apertures 

We use Monte Carlo simulations of the observed stellar ve- 
locities and corresponding errors to investigate the recov- 
ery of their average kinematics. We mimic the stellar ve- 
locity observations by randomly drawing from an assumed 
intrinsic Gaussian velocity distribution, with given mean 
velocity Vq and velocity dispersion uq. This set of intrin- 
sic velocities, is then 'instrumentally' broadened by adding 
to each velocity a random drawing from a Gaussian with 
zero mean and the velocity error as standard deviation. 
These velocity errors are simulated by randomly drawing 
from the observed velocity error distribution (right panels 
of Figure |5J) . For the latter we use the rejection method 
(§ 7.3 of Press et al. 1992), with a Lorentzian distribu- 
tion as comparison function. In this way, we create, for a 
given number of stars, 500 sets of simulated velocities and 
corresponding errors. 

Next, we use the maximum likelihood method of 
Appendix E] to calculate the mean velocity and velocity 
dispersion for each simulated set separately. In Figure lA.ll 
we compare the (biweight 10 ) mean (filled circles) of these 
500 mean velocity and velocity dispersion measurements 
with Vq and oq (horizontal lines) of the given intrinsic 
Gaussian velocity distribution. The error bars are the 
(biweight) standard deviation of the kinematic measure- 
ments, and indicate the precision with which the kine- 
matics can be measured, given the observed velocity er- 
ror distribution. The precision increases with increasing 
number of stars per bin. The precision also increases with 
decreasing intrinsic mean velocity dispersion ctq. To re- 
move the latter dependency, we give relative kinematic 

10 The biweight mean and biweight standard deviation (e.g. 
Andrews et al. 1972; Beers et al. 1999) are robust estimators 
for a broad range of non- Gaussian underlying populations and 
are less sensitive to outliers than other moment estimators. 
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Fig. A.l. Recovery of maximum-likelihood-estimated 
kinematics from proper motions (top panel) and line-of- 
sight velocities (bottom panel). For a given number of 
stars per aperture, velocities and corresponding errors are 
simulated by randomly drawing from an intrinsic Gaussian 
distribution with mean velocity Vo and velocity disper- 
sion (To, broadened by velocity errors randomly drawn 
from the observed velocity error distributions (left panels 
FigureEJ. Each filled circle with error bar shows the mean 
and standard deviation of the measured kinematics from 
500 such simulations. As a compromise between lower pre- 
cision (larger error bars) for a small number of stars per 
aperture, and lower spatial resolution (larger bins) for a 
larger number of stars, we choose to have between 50 and 
100 stars per bin. 



measurements and corresponding errors, i.e., divided by 
the (arbitrarily) chosen value for o~q. 

Both the mean velocity and velocity dispersion are re- 
covered well. To obtain a better precision, we can increase 
the number of stars per aperture, but at the same time the 
spatial resolution decreases, as we have to increase the size 
of the apertures. We find that between 50 and 100 stars 
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per aperture is a good compromise. For the proper motions 
this implies a (relative) precision for the mean velocity V 
and velocity dispersion a of respectively AV/a ~ 0.12 and 
Aa/a ~ 0.09. For the line-of-sight velocities we find simi- 
lar values, respectively AV/a ~ 0.12 and Aa/a ~ 0.08. 

Given the average proper motion dispersion of about 
0.5 masyr -1 for w Cen (§ I7.2|l . this means we expect to 
measure the mean proper motion and dispersions with an 
average (absolute) precision of respectively 0.06 masyr -1 
and 0.05 masyr -1 . Similarly, with an average line-of-sight 
velocity dispersion of about 14 kms -1 for to Cen, we ex- 
pect to measure the mean line-of-sight velocity and disper- 
sion with an average precision of respectively 1.7 kms -1 
and 1.1 kms -1 . 

Indeed, the average of the uncertainties in the kine- 
matics given in Table |3] and 21 are consistent with these 
expectations. Moreover, as predicted, the decrease in the 
uncertainties with radius is proportional to the decrease in 
dispersion. In other words, if we divide the uncertainties 
by the corresponding dispersions, we find nearly constant 
(relative) precisions, AV/a ~ 0.11 and Aa/a ~ 0.08 for 
both proper motions and line-of-sight velocities, consistent 
with the above simulated precisions. 

To enhance the signal-to-noise of the observations, we 
first reflect all measurements back to the first quadrant 
(x' > 0, y' > 0). We exploit the fact that for an axisym- 
metric object, the proper motions in the x'-direction are 
symmetric in the projected minor axis, while the proper 
motions in the j/'-direction as well as the line-of-sight ve- 
locities are symmetric in the projected major axis. Since 
our models are intrinsically axisymmetric, it is equivalent 
to fit either to the original or to the reflected data. 

We use a polar grid of apertures on (the first quadrant 
of) the plane of the sky to better approximate the shape of 
photometric and kinematic observations. Every aperture 
is characterised by its central radius r$ > and angle 
0° < 6*o < 90°, together with its radial and angular width, 
denoted by Ar and AO, respectively. We construct the 
polar grids such that each aperture has (at least) 50 stars, 
together with the requirement that the apertures are as 
'round' as possible in the sense that Ar w roA6>. The 
latter avoids (very) radial or angular elongated apertures, 
which would include stars from (very) different positions, 
with probably different (kinematical) properties than the 
stars near the centre of the aperture. 

Appendix C: Simple distance estimate 

The most straightforward way to obtain a dynamical dis- 
tance estimate is from the ratio of the line-of-sight veloc- 
ity dispersion <7i os and the proper motion velocity disper- 
sion a pm for spherically symmetric objects (e.g. Binney & 
Tremaine 1987, p. 280) 



D (kpc) 



trios (kms : ) 



(C.l) 



718 stars for which all three velocity components are mea- 
sured, we find for the two mean proper motion disper- 
sion components a x > — 0.58 ± 0.02 masyr -1 and a y > = 
0.55 ±0.02 masyr -1 , and for the mean line-of-sight veloc- 
ity dispersion a z > — 12.3 ± 0.3 kms -1 . Substituting the 
latter value together with the average proper motion dis- 
persion in (|C.lfl . we obtain a distance of D = 4.6±0.2 kpc. 

This value is below the canonical distance D = 5.0±0.2 
(Harris et al. 1996). The above simple distance estimate 
is not valid for to Cen, which is not spherically symmetric. 
Moreover, although the above average values for a x i and 
ay are just consistent with each other, from the left panel 
of Figurc lB~T1 it is clear that the profile of the mean proper 
motion dispersion profile of a x * (green) lies systematically 
above that of oy (red). A non-spherical anisotropic model 
is needed to explain these observations. Here we consider 
a simple model with constant anisotropy. 

If we make the (ad-hoc) assumption that the ve- 
locity ellipsoid is oblate with intrinsic semi-axis lengths 
a x = a y = a and a z = q vc a (all in kms -1 ), where q vo 
is the average intrinsic flattening, the observed velocity 
dispersions are given by 



fx' 

ov 



a / masyr 1 , 

q vc a / 4.74D masyr -1 , 

[l - {1 - q' v 2 c ) cot 2 i] 1/2 a kms -1 , 



(C.2) 



4.74 er pm (masyr x ) ' 

Using, from the 2295 selected stars with proper motions 
and 2163 selected stars with line-of-sight velocities, the 



where we have used eq. (J2J and the relation q 2 sin 2 i — 
q' 2 — cos 2 i. Using the best-fit value for Z?tani of 5.6 
kpc (§ I4.5JI . we eliminate the inclination i. Next, by fit- 
ting the ratios of the line-of-sight velocity dispersion over 
the proper motion dispersion components, ov/ov and 
0V / cry' , to the observations in the left panel of Figure 
we determine the best-fit values for the remaining two free 
parameters: the distance D and the (projected) flattening 
of the velocity ellipsoid q' vc . 

Since we use the full dispersion profiles and we allow for 
an anisotropic velocity distribution, this simple way to ob- 
tain a dynamical distance estimate goes beyond the above 
spherical symmetric approach. If q' vc = 1 in eq. ljC.2fl . we 
recover this spherical symmetric approach in which both 
ratios are equal and the distance follows from eq. ilC.l|) . 

We show in the right panel of Figure IB. II the A\ 2 
contours for a range of q' vc and D. The overall minimum, 
indicated by a cross, corresponds to the best-fit values 
q' ve = 0.92 ± 0.05 and D = 4.54 ± 0.14 kpc. The isotropic 
case (q' ve = 1) is excluded at about the 95.4%-level. The 
best-fit (projected) flattening of the velocity ellipsoid is 
less than the average observed flattening q 1 — 0.879±0.007 
(hashed region) from the stellar photometry of to Cen 
(Geyer et al. 1983), although an equivalent value is not 
excluded (at the 68.3%-level). The velocity distribution is 
expected to be less flattened, since it traces more directly 
the potential, which in general is rounder than the light 
distribution (see e.g. p. 48 of Binney & Tremaine 1987). 

If we only fit the ratio ov /ov, the green dashed curve 
shows the best-fit distance at given flattening. While in 
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this case the distance increases with flattening, almost ex- 
actly the opposite happens if we only fit the ratio <j z i /ay 
(red dashed curve). Simultaneously fitting both ratios does 
not provide a good fit (the x 2 value is significantly larger 
than the number of degrees of the freedom) and the result- 
ing best-fit distance (black dashed curve) of about 4.5 kpc 
is significantly below the canonical distance of 5.0 kpc. 

We conclude that both the simple distance estimate 
l)C.l(l and the above constant-anisotropy axisymmetric 
model are not valid for u) Cen and underestimate its 
distance. To explain the observed kinematics of u> Cen 
and obtain a reliable distance estimate, one needs a non- 
spherical dynamical model with varying anisotropy, like 
the Schwarzschild modelling technique used in this paper. 

Finally, it is interesting to note that Peterson & 
Cudworth (1994) reported rotation in the line-of-sight ve- 
locities and proper motions of M22, and found that its dy- 
namical distance increased slightly after an approximate 
correction based on a comparison of dispersion profiles. 
Peterson et al. (1995) saw no evidence for rotation in M4, 
but did note that their resulting dynamical distance was 
smaller than the canonical value. Both studies used the 
simple distance estimates described here. It will be inter- 
esting to reanalyse these clusters with the comprehensive 
method we have presented in this paper. 
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Fig. 20. Kinematics of different components in the distribution function of our best-fit model for w Cen. From left to 
right: full distribution function, main inner component, main outer component and separate disk component between 
1 and 3 arcmin (see text for details). From top to bottom: spatial distribution, mean velocity fields in the direction 
of the major x'-axis, the minor y'-axis and the line-of-sight z'-axis, and mean velocity dispersion profiles. The radial 
dispersion ctr< (green) and tangential dispersion ag> (red) are on the plane of the sky and ov (blue) is the line-of-sight 
dispersion. 
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Fig. B.l. Left panel: velocity dispersion profiles calculated along concentric rings. Assuming the canonical distance 
of 5 kpc, the profiles of the proper motion components in the aZ-direction (green) and y'-direction (red) are converted 
into the same units of kms -1 as the line-of-sight profile in the z'-direction (blue). The black horizontal lines indicate 
the corresponding scale in masyr . Below the profiles, the mean velocity error per ring is indicated. Right panel: 
Assuming an oblate velocity ellipsoid with constant (projected) flattening, the ratio of the line-of-sight over the proper 
motion velocity dispersion profiles yields an estimate for the dynamical distance D. The best-fit values correspond 
to the minimum (cross) in the A% 2 contour plot, where the inner three contours are drawn at the 68.3%, 95.4% and 
99.7% (thick contour) levels, and subsequent contours correspond to a factor of two increase in A% 2 . For increasing 
flattening of the velocity ellipsoid, starting with the isotropic case on the left axis, the green (red) dashed curve shows 
the corresponding best-fit distance if only the profile of the proper motion in the x'-direction (j/'-direction) is used, 
and the black dashed curve if both are used. The observed flattening from the stellar photometry (Geyer et al. 1983) 
is indicated by the hashed region. 



